Welcome to the Subcellular workflow documentation!
(Under construction - last updated 2021/11/03)
This workflow has been developed to tackle the challenge of building and analyzing biochemical pathway models, combining pre-existing tools and custom-made software. (Santos et al. 2020) (Preprint)
At the root of our implementation is the SBtab format (Lubitz et al. 2016), a file that can store biochemical models and associated data in an easily readable and expandable way.
We have also developed tools to convert the SBtab format into several formats that can be used in MATLAB®, NEURON, STEPS and COPASI.
Using MATLAB® we have developed custom scripts for parameter estimation, and global sensitivities analysis, as well as diagnostics tools that can be used for model development. The global sensitivity analysis algorithm is modified from Halnes et al 2009.
We demonstrate all these features using three example models, the main one being a modified version of the D1 MSN subcellular cascade model from Nair et al. 2016.
Code and files to run these models in different simulators:
- For Neuron
- For the Subcellular application (STEPS)
- Copasi
Features:
Wrapper for model simulation in MATLAB® (Matlab/Run_main.m)
Analysis of selected parameter sets, using MATLAB® (Matlab/Run_main.m)
Parameter optimization, using MATLAB® (Matlab/Run_main.m)
Global Sensitivity analysis, using MATLAB® (Matlab/Run_main.m)
Conversion tools:
- SBtab (.xlsx,.xls) to SBtab (.tsv), using MATLAB® (Matlab/Run_main.m)
- SBtab (.xlsx) to MATLAB® SimBiology® (.m, .sbproj), using MATLAB® (Matlab/Run_main.m)
- MATLAB® SimBiology® to SBML (.xml), using MATLAB® (Matlab/Run_main.m)Needs to be fixed with our R script (https://github.com/a-kramer/simbiology-sbml-fix)
- SBtab (.tsv) to VFGEN (.vf), using R (https://github.com/a-kramer/SBtabVFGEN)
- SBtab (.tsv) to Mod (.mod), using R (https://github.com/a-kramer/SBtabVFGEN)
- SBtab (.tsv) to SBML (.xml), using R (https://github.com/a-kramer/SBtabVFGEN)
Sections of the workflow in external repositories
Conversion tools:
Implemented models
Compatibility
Subcellular workflow MATLAB® code is compatible with MATLAB® 2020a or above running on Microsoft Windows, macOS and Linux.
Matlab® packages needed:
Optimization Toolbox™
Statistics and Machine Learning Toolbox™
Fuzzy Logic Toolbox™
Financial Toolbox™
Global Optimization Toolbox
SimBiology®
Parallel Computing Toolbox™
SBtab
This is an overview of the SBtab syntax that is used in our workflow. This contains a list of the sheet names and subfields that are read by our software. The second column in the first row of each sheet must include “TableName=’sheet name’”. Fields that are not mentioned here are not used in the latest workflow and are not imported but might be added to future releases. Additional information on how to use SBtab can be found in https://www.sbtab.net/sbtab/default/documentation.html.
Defaults
!ID - Identifier code for the entries, it should consist of the letter “D” followed by an integer, should start at 1.
!Name - Name of the default variable being defined, we have used time, volume, substance, lenght, and area in our files.
!Unit - Unit of the default variable, eg. time - second.
Compartment
!ID - Identifier code for the entries, it should consist of the letter “V” followed by an integer, should start at 1.
!Name - Name of the compartment.
!Unit - Compartment volume unit, usually liter.
!Size - Size of the compartment in units defined in the unit column.
Compound
!ID - Identifier code for the entries, it should consist of the letter “S” followed by an integer, should start at 0.
!Name - Name of the compound. We advise to use recognizable compound names and separate complexes of multiple compounds with ‘_’ (e.g. A_B). Must start with a letter.
!Unit - Unit for the compound, usually nanomole or nanomole/liter. If it is not a default MATLAB |Reg| unit it should be added to it before trying to run the scripts.
!InitialValue - Default initial values for the compounds before the equilibration step. These are usually overriden by the values Si in the experiments sheet
!IsConstant - assigns a binary value, either TRUE or FALSE, depending on whether the value of a particular compound should stay constant throughout the simulations. Note that the input compound should remain constant.
!Assignment - assigns a binary value, either TRUE or FALSE.
!Interpolation -
!Type - Type of the reaction, e.g. “kinetic”.
!Location - Compartment in which a particular compound is present.
Reaction
!ID - Identifier code for the entries, it should consist of the letter “R” followed by an integer, should start at 0.
!Name - Reaction name. Must start with a letter and it is advisable to include the reaction index, e.g. ReactionFlux0.
!KineticLaw - Reaction kinectic law. Needs to include the precise mass action kinetic law with compound and parameter names from corresponding SBtab table. For a reaction ‘A + B <=> A_B’ with a forward reaction rate of kf_R0 and backward reaction rate of kr_R0 the formula would look like ‘kf_R0*A*B-kr_R0*A_B’.
!IsReversible - Bollean identifying the reversibility of the reaction, it is TRUE for reversible and FALSE for irreversible reactions.
!Location - Compartment in which a particular reaction is taking place.
!ReactionFormula - Chemical formula of the reaction, should be written in the form ‘A + B <= > A_B’.
Parameter
!ID - Identifier code for the entries, it should consist of the letter “K” followed by an integer, should start at 0.
!Name - Name of the parameter. We followed the convention of using ‘kf’ for forward reactions rates and ‘kr’ for reverese rates, followed by the reaction !ID, e.g. ‘kf_R0’. These names can be arbitrary but they need to coincide with whatever is defined in the Reaction kinectic law
!Unit - The units of the parameter.
!DefaultValue - Parameter value in linear space.
!Value:linspace - Parameter value in linear space.
!Value:log2 - Parameter value in log base 2.
!Value:log10 - Parameter value in log base 10.
!Location - Compartment in which a particular parameter is governing a reaction.
!Comment - Could be any plain text , we used it as a handy way of determining which reaction the parameter is involved in. We advise to use the following syntax ‘kf_AXB__A_B’ and ‘kr_AXB__A_B’ for respectively forward and backward rates of a reaction ‘A + B <= > A_B’.
Expression
Compounds which are defined by expressions.
!ID - Identifier code for the entries, it should consist of the letters “Ex” followed by an integer, should start at 0.
!Name - Name of the compound defined by the expresion.
!Formula - Formula assigned to the compound, it should use the compound names used in the compound sheet and, if needed, constant names from the constant sheet.
!Unit - Concentration unit for the Compound.
!Location - Compartment in which a particular compound is located.
Output
Compounds used as outputs in experimental data.
!ID - Identifier code for the entries, it should consist of the letter “Y” followed by an integer, should start at 0.
!Name - Name used to identify the output compound, when an existing compound needs to be measured we usually use “compound_name”_out.
!ErrorName - Name of the error of the output compound. It should start with ‘SD’ (referring to standard deviation) followed by the output ID, e.g. SD_Y0.
!ErrorType - Type of error for the output compound, we have used ‘abs+rel random noise (std)’.
!Unit - Concentration unit for the output compound.
!ProbDist - Probability distribution type of the measured output in an experimental setting, e.g. ‘normal’.
!Location - Compartment in which a particular output compound is located.
!Formula - Formula that links the experimental measured output to the compounds in the model. Usually the experimental measurement corresponds to a sum of compounds existing in the model but ratios are also common.
Experiments
Each column corresponds to one experiment for which there is a separate sheet.
!ID - Identifier code for the entries, it should consist of the letter “E” followed by an integer, should start at 0.
!Name - Name used to identify the experiment, we advise using the the word ‘Experiment’ followed by the experiment index.
>Output - Should list all the output !ID’s, i.e. Y’s followed by their indices and separated by commas.
>Si- List of the various compounds of the model that have starting amounts other than 0, the same !ID as in the coumpound table should be used.
In a model with 2 experiments and 5 compounds A,B,C,D,E with IDs S0,S1,S2,S3,S4 respectively
Experiment1 with compounds starting amounts A=0,B=1,C=2,D=0,E=3
Experiment2 with compounds starting amounts A=1,B=0,C=1,D=0,E=4
Four entries should be included as exemplified bellow:
>S0
>S1
>S2
>S4
Experiment1
0
1
2
3
Experiment2
1
0
1
4
Note that D/S3 is omitted because the starting value is 0 in all experiments.
!SimTime - Simulation time for a particular experiment.
!Normalize - Normalizations to be performed to the outputs are defined here.
Ei
Corresponds to individual experiments. It should have the name of the experiment IDs used in the experiments sheet.
!ID - Identifier code for the entries, it should consist of the letters “EiT” followed by an integer, should start at 0.
!Time - Time series data, this should include a list of all the time points during which the corresponding output data points were sampled.
`>Yi- Compound to be measured, corresponds to the !ID of the output sheet. It should have the experimental (or simulated from another model) data.
`SD_Yi- Error of the compound to be measured, corresponds to the !ErrorName of the output sheet. It should have the experimental (or simulated from another model) data.
EiI
Corresponds to individual experiments. It should have the name of the experiment IDs used in the experiments sheet, with an i in between E and the experiment number.
!ID - Identifier code for the entries, it should consist of the letters “EiIT” followed by an integer, should start at 0.
!Input_Time_Si- Time series of the inputs to the model, this should include a list of all the time points during which the corresponding input data points were sampled. To produce simple step inputs, only the time points during which a change in concentration is happening can be included. To produce more complicated input curves, more time points are needed to represent the shape of the curve.
`>Si- Compound that is being changed as input to the model, corresponds to the !ID in the compound table. This column should represent the sampled concentration data points corresponding to each time point.
MATLAB®
The MATLAB® section of this workflow has been developed to facilitate model building and rapid iteration between different versions of a model. In this workflow we use one main script, “Run_main.m”, that calls all the relevant functions to be used. To run the MATLAB® code the Subcellular Workflow repository should be added to the MATLAB® path. Running the script “Run_main.m” generates prompts in the MATLAB® terminal window with a request to the user, to choose between a number of different options. These prompts allow the user to choose:
The model to use (from all the models that are in the “Matlab/model” folder).
Note that the very first time you run a model you have to add the folder for that specific model from its home repository into the “Matlab/model” folder(e.g. copy the folder “Model_nair_2016” from its repository to “Matlab/model/” ).For implemented models so far go to the following links:The analysis to be performed, with the following options:
Reproduction of a previous analysis
This option can be used to re-do an analysis that has previously been performed. This is useful for reproducibility and in the case of the code getting updated with extra funcionalities. The user should specify the analysis file that they want to use, examples are provided in the each model repository.
Reproduction of the plots of a previous analyis
Similar to the previous option but here only the plots are re-done.
Import model files
Creation of the model files and folder that are needed to run the model in Matlab, the creation of a folder with the model in .tsv format (one tsv file for each excel sheet of the original SBtab), as well as the conversion of the model to the SBML format (.xml).
Examples of the output recieved when the different models are run through the workflow can be found on the respective model repository in the directory “Matlab/Results/Results/Examples”, or by following these links:
In order to gain a better understanding on how the code works, there are detailed pages for the following:
Scripts - The script that we use;
Functions - All the custom functions we have built; this is directed to anyone that wants to develop or iterate the code;
Settings file - The master configuration file, where we describe everything that can be modified by the user without changing any code;
Results - Explanation of all the files containing relevant results that are generated after running the built-in analysis of the code;
Model files and folders - Description of all the files and folders that are generated when the model is imported from SBtab into relevant files, used by the rest of the MATLAB® code.
Diagnostics
The specifications of the diagnostics operations require user input in the settings file. The toolkit imports the model stored in SBtab format (with the name specified by stg.sbtab_excel_name variable from a folder chosen by the user terminal prompts) to a .sbproj file, and saves it in the subfolder called Data along with the imported data and inputs in a .mat format. Once the model has been imported, the import function can be disabled in the import section of the settings file for further procedures.
The parameter sets that are specified in the diagnostics section of the settings file are then used in model simulations with the input specifications (e.g. time series data for relevant input species) in the Experiment Input (EI) tables in the SBtab file, calculate the error between the simulation results and the experimental data, and plot the error scores as well as the comparative traces of the simulation results and the experimental data. At least one parameter set is currently required in the settings file. Experiments of interest can be specified by the stg.exprun variable in the analysis section of the settings file.
Parameter Estimation
Parameter estimation is performed by various MATLAB optimization algorithms. The number of parameters to estimate, possible thermodynamic constraints (can be determined by a standalone script), and upper and lower bounds can be specified in the model section of the settings file. In addition, the parameter indices and the best available parameter sets can specified in the diagnostics section. Optimization algorithm and optimization settings can be found in the optimziation section of the settings file and simulation settings in the simulation section.
Global Sensitivity Analysis
Global Sensitivity Analysis can be performed when a user is interested in parameter distribution rather than single point values, and how sensitive a specific output is towards perturbations in different parameters. Instructions on what settings are required to be specified can be found in the Global Sensitivity Analysis section of the settings file.
Scripts
This is the entry point for our code, it calls all other relevant functions.
Run_main
Code
1% Script to import sbtab and run the analyis 2% clear functions 3 4%Get the date and time 5date_stamp = string(year(datetime)) + "_" + ... 6 string(month(datetime)) + "_" + string(day(datetime))... 7 + "__" + string(hour(datetime)) + "_" + string(minute(datetime))... 8 + "_" + string(round(second(datetime))); 9 10% Get the folder where the MATLAB code and models are located 11Matlab_main_folder = fileparts(mfilename('fullpath'))+"/"; 12Matlab_main_folder = strrep(Matlab_main_folder,"\","/"); 13addpath(genpath(Matlab_main_folder)); 14mmf.main = Matlab_main_folder; 15 16% Name of the various analysis that can be run with this code 17analysis_options = ["Diagnostics","Parameter Estimation",... 18 "Global Sensitivity Analysis","PLA","Reproduce a previous analysis",... 19 "Reproduce the plots of a previous analysis","Import model files"]; 20 21% Code for choosing the model and loading the settings files 22[stg,rst,sb] = f_user_input(mmf,analysis_options); 23 24% Get the folder structure used for the model files 25[mmf] = default_folders(stg,mmf,date_stamp); 26 27% Runs the import scripts if chosen in settings 28if stg.import 29 [stg,sb] = f_import(stg,mmf); 30else 31 % Creates a struct based on the sbtab that is used elswhere in the code 32 % and also adds the number of experiments and outputs to the settings 33 % variable 34 if isempty(sb)% check needed for plot reproduction 35 [stg,sb] = f_generate_sbtab_struct(stg,mmf); 36 end 37end 38 39% Runs the Analysis chosen in settings 40if any(contains(analysis_options(1:4),stg.analysis)) 41 rst = f_analysis(stg,stg.analysis,mmf,analysis_options); 42end 43% Save Analysis results if chosen in settings 44if stg.save_results 45 f_save_analysis(stg,sb,rst,mmf) 46end 47 48% Plots the results of the analysis, this can be done independently after 49% loading the results of a previously run analysis 50if stg.plot 51 f_plot(rst,stg,mmf) 52 % Save plots results if chosen in settings 53 if stg.save_results 54 f_save_plots(mmf) 55 end 56end
This is the main script from the MATLAB® portion of the workflow. Depending on the configurations on the settings file and choices on the user facing prompts it can call functions to:
Perform conversions of the SBtab:
SBtab( .xlsx) to SBtab (.tsv)
SBtab (.xlsx) to MATLAB® SimBiology® (.m, .sbproj)
MATLAB® SimBiology® to SBML (.xml)
Perform analysis on the model:
Diagnostics
Parameter Estimation
Global Sensitivity Analysis
It can also reproduce a the calculations of a previous analysis or just its plots.
Functions
The MATLAB® functions used in this workflow are divided acording to their role, we have:
Functions to import the model to MATLAB® and generate model specific files and functions.
Simulation and scoring functions
Functions relating to the simulation and scoring of the model.
Functions for the analysis that we can perform on the model.
Functions to plot the result of the different analyses.
General purpose functions that are usually used by other functions.
Setup and Import
f_user_input
Code
1function [stg,rst,sb] = f_user_input(mmf,analysis_options) 2 3persistent last_SBtab_date 4persistent last_model_folder 5persistent last_settings_file_text 6persistent last_settings_file_date 7persistent last_analysis_text 8 9Matlab_main_folder = mmf.main; 10 11rst = []; 12sb = []; 13functions_cleared = false; 14 15% Get the folder of the model 16model_folder_general = Matlab_main_folder + "Model/"; 17last_choice = last_model_folder; 18prompt = "What model folder should be used?\n"; 19[model_folder,last_model_folder] =... 20 choose_options(model_folder_general,prompt,last_choice); 21 22model_name_specific = string(model_folder); 23folder_model_specific = model_folder_general + "/" + model_name_specific; 24 25prompt = "\nWhat analysis should be performed?\n"; 26last_choice = last_analysis_text; 27[analysis_text,last_analysis_text] =... 28 parse_choices(prompt,analysis_options,last_choice); 29 30% analysis_n = find(contains(analysis_options,analysis_text)); 31 32% Check if an analysis was chosen 33if any(contains(analysis_options([1:4,7]),... 34 analysis_text)) 35 36 %Get the Setting file to be used 37 settings_folder = folder_model_specific + "/Matlab/Settings"; 38 last_choice = last_settings_file_text; 39 prompt = "\nWhat file should be used as settings?\n"; 40 [settings_file_text,last_settings_file_text] =... 41 choose_options(settings_folder,prompt,last_choice); 42 43 settings_file = strrep(settings_file_text,".m",""); 44 45 % Add the default settings to the struct 46 stg_add_default = eval("default_settings()"); 47 48 f = fieldnames(stg_add_default); 49 for i = 1:length(f) 50 stg.(f{i}) = stg_add_default.(f{i}); 51 end 52 53 % Add chosen settings to the struct overwriting defaults when 54 % appropriate 55 [stg_add] = eval(settings_file + "()"); 56 57 f = fieldnames(stg_add); 58 for i = 1:length(f) 59 stg.(f{i}) = stg_add.(f{i}); 60 end 61 62 % Check if the date of the settings file changed, if so clear functions 63 listing = dir(settings_folder); 64 for n = 1:size(listing,1) 65 if matches(settings_file_text,listing(n).name,"IgnoreCase",true) 66 settings_file_date = listing(n).date; 67 end 68 end 69 70 [last_settings_file_date,functions_cleared] =... 71 compare_last(settings_file_date,last_settings_file_date,... 72 functions_cleared); 73 74 % Check if the name of the settings file changed, if so clear functions 75 [last_settings_file_text,functions_cleared] =... 76 compare_last(settings_file_text,last_settings_file_text,... 77 functions_cleared); 78 79 % Check if the date of the SBtab changed, if so clear functions 80 listing = dir(folder_model_specific); 81 82 for n = 1:size(listing,1) 83 if matches(stg.sbtab_excel_name,listing(n).name,"IgnoreCase",true) 84 sbtab_date = listing(n).date; 85 end 86 end 87 [last_SBtab_date,~] =... 88 compare_last(sbtab_date,last_SBtab_date,functions_cleared); 89 90 % Store the name of the chosen analysis in the settings struct 91 stg.analysis = analysis_text; 92 93 if contains(analysis_options(7),analysis_text) 94 stg.import = true; 95 stg.save_results = false; 96 stg.plot = false; 97 end 98elseif any(contains(analysis_options(5:6),analysis_text)) 99 100 % Get the folder of the Analysis that should be reproduced 101 folder_results = folder_model_specific + "/Matlab/Results"; 102 103 last_choice = []; 104 prompt = "\nWhat analysis should be reproduced?\n"; 105 106 [r_analysis_text,~] =... 107 choose_options(folder_results,prompt,last_choice); 108 109 folder_results_specific = folder_results + "/" + ... 110 r_analysis_text; 111 112 last_choice = []; 113 prompt = "\nWhen was this analysis run originaly?\n"; 114 115 [r_analysis_date_text,~] =... 116 choose_options(folder_results_specific,prompt,last_choice); 117 118 folder_results_specific_date = folder_results_specific + "/" + ... 119 r_analysis_date_text; 120 121 % Load the settings file and the SBtab struct 122 load(folder_results_specific_date + "/Analysis.mat","stg","sb") 123 124 % Set inport to false since we don't want to overwrite anything 125 stg.import = false; 126 127 % If the reproduction of an analysis is chosen clear the functions 128 % because the settings most likely changed 129 if contains(analysis_options(5),analysis_text) 130 f_functions_to_clear() 131 end 132 133 % I the reproduction of the plots of an analyis is chosen make sure 134 % we tell the code to produce plots and also load the results that were 135 % previously obtained 136 if contains(analysis_options(6),analysis_text) 137 stg.plot = true; 138 load(folder_results_specific_date + "/Analysis.mat","rst") 139 end 140end 141 142% Set the chosen model folder in the settings struct 143stg.folder_model = model_name_specific; 144end 145 146function [choice,last_choice] = choose_options(folder,prompt,last_choice) 147 148listing = dir(folder); 149 150for n = size(listing,1):-1:1 151 if any(matches(listing(n).name,[".","..","Place models here.txt"])) 152 listing(n)= []; 153 end 154end 155 156for n = 1:size(listing,1) 157 options(n) = string(listing(n).name); 158end 159 160[choice,last_choice] = parse_choices(prompt,options,last_choice); 161end 162 163function [choice,last_choice] = parse_choices(prompt,options,last_choice) 164 165for n = 1:size(options,2) 166 prompt = prompt + "\n" + n + ": " + options(n); 167end 168 169if ~isempty(last_choice) 170 if any(contains(options,last_choice)) 171 prompt = prompt + "\n\nPress enter to use " + last_choice; 172 else 173 last_choice = []; 174 end 175end 176 177prompt = prompt + "\n"; 178 179i = input(prompt); 180 181if isempty(i) 182 choice = []; 183elseif i > 0 && i < size(options,2)+1 184 choice = options(i); 185 disp("The option chosen was: " + choice) 186else 187 prompt = "Please choose from the provided options"; 188 [choice,last_choice] = parse_choices(prompt,options,last_choice); 189end 190 191if isempty(choice) 192 if ~isempty(last_choice) 193 choice = last_choice; 194 disp("The option chosen was: " + last_choice) 195 else 196 prompt = "Please choose from the provided options"; 197 [choice,last_choice] = parse_choices(prompt,options,last_choice); 198 end 199else 200 last_choice = choice; 201end 202end 203 204function [previous,f_cleared] = compare_last(current,previous,f_cleared) 205 206if ~isempty(previous) 207 if ~contains(current,previous) 208 if f_cleared == false 209 disp("Settings file changed, clearing functions") 210 f_functions_to_clear() 211 f_cleared = true; 212 end 213 end 214end 215previous = current; 216end
It prompts the user to choose the model to run, the settings file to use, and the Analysis to perform.
f_import
Code
1function [stg,sb] = f_import(stg,mmf) 2 3Model_folder = mmf.model.main; 4 5disp("Generating model files and folder from SBtab") 6 7% Create needed folders 8[~,~] = mkdir(mmf.model.data.main); 9[~,~] = mkdir(mmf.model.input_functions.main); 10[~,~] = mkdir(mmf.model.tsv.model_name); 11[~,~] = mkdir(mmf.model.data.model_exp.main); 12 13% Creates a .mat and a tsvs from the SBtab file 14f_excel_sbtab_importer(mmf); 15 16addpath(genpath(Model_folder)); 17 18% Creates a struct based on the SBtab that is used elswhere in the code and 19% also adds the number of experiments and outputs to the settings struct 20[stg,sb] = f_generate_sbtab_struct(stg,mmf); 21 22%Create the model and input output structure from sbtab.mat 23 24% Saves the model in .mat, .sbproj and .xml format, while also creating a 25% file whith the data to run the model in all different experimental 26% settings defined in the SBtab 27f_sbtab_to_model(stg,sb,mmf) 28 29% Creates code that loads the inputs of each experiment into a .mat file, 30% and creates the code to read this inputs at runtime when the experiments 31% are being simulated, all this generated code is stored on the Input 32% functions folder 33f_setup_input(stg,mmf) 34 35%Creates three .mat files for each experiment, with all the added rules, 36%species and parameters needed depending on the inputs and outputs 37%specified on the SBtab, one for the equilibrium simulation run, one for 38%the deffault run, and one for a more detailed run. 39f_build_model_exp(stg,sb,mmf) 40disp("Model files and folders generated successfully") 41end
Creates the necessary folders inside the model folder. Calls subfunctions that convert the SBtab from an Excel into MATLAB® files useful for the workflow, TSVs and a SBML.
f_excel_sbtab_importer
Code
1function f_excel_sbtab_importer(mmf) 2 3Source_sbtab = mmf.model.sbtab; 4Matlab_sbtab = mmf.model.data.sbtab; 5 6% Get the total number of sheets in the SBTAB 7sheets = sheetnames(Source_sbtab); 8 9% Try to run the import the sheets in multicore, depending on the version 10% of excel this migth not work 11try 12 parfor i = 1:size(sheets,1) 13 sbtab_excel{i} = impexp (i,mmf); 14 end 15catch 16 for i = 1:size(sheets,1) 17 sbtab_excel{i} = impexp (i,mmf); 18 end 19end 20 21% Save the SBTAB tables in .mat format 22save(Matlab_sbtab,'sbtab_excel'); 23disp("SBtab with " + size(sheets,1) + " sheets parsed successfully") 24end 25 26function sbtab_excel = impexp (i,mmf) 27 28Source_sbtab = mmf.model.sbtab; 29tsv_name_folder = mmf.model.tsv.model_name; 30 31% Import the SBTAB to a cell with a sheet per cell 32sbtab_excel = readcell(Source_sbtab,'sheet',i); 33 34% Replace "ismissing" values with empty spaces 35mask = cellfun(@ismissing, sbtab_excel,'UniformOutput',false); 36mask = cellfun(@min, mask); 37mask = logical(mask); 38sbtab_excel(mask) = {[]}; 39 40% Get name for tsv that is going to be exported 41field = regexp(sbtab_excel{1,2},"TableName='[^']*'",'match'); 42field = string(replace(field,["TableName='","'"," "],["","","_"])); 43 44%Export the tsv 45cell_write_tsv(tsv_name_folder + field + ".tsv",sbtab_excel) 46end 47 48function cell_write_tsv(filename,origCell) 49 50% save a new version of the cell for reference 51modCell = origCell; 52% assume some cells are numeric, in which case set to char 53iNum = cellfun(@isnumeric,origCell); 54 55% Replace numeric sells with cell strings 56for n = 1:size(iNum,1) 57 for m = 1:size(iNum,2) 58 modCell(n,m) = cellstr(num2str(origCell{n,m})); 59 end 60end 61 62%Save the file that only as strings in each cell 63modCell = transpose(modCell); 64 65[rNum,cNum] = size(origCell); 66frmt = repmat([repmat('%s\t',1,cNum-1),'%s\n'],1,rNum); 67fid = fopen(filename,'wt'); 68fprintf(fid,frmt,modCell{:}); 69fclose(fid); 70end
Loads the information in the SBtab and creates a . mat file that contains the sbtab and TSVs corresponding to all the SBtab tabs.
Inputs - stg
Saves
.mat file containing the SBtab in the “Model/Data” folder
TSVs containing the SBtab in the “Model/tsv” folder
f_generate_sbtab_struct
Code
1function [stg,sb] = f_generate_sbtab_struct(stg,mmf) 2 3Matlab_sbtab = mmf.model.data.sbtab; 4 5if isfile(Matlab_sbtab) 6 7 load(Matlab_sbtab,'sbtab_excel'); 8 9 sb = f_get_sbtab_fields(sbtab_excel); 10 11 stg.expn = size(sb.Experiments.ID,1); 12 stg.outn = size(sb.Output.ID,1); 13end 14end 15 16function sb = f_get_sbtab_fields(sbtab_excel) 17for n = 1:size(sbtab_excel,2) 18 19 if ~isempty(sbtab_excel{1,n}{1,2}) 20 21 field = regexp(sbtab_excel{1,n}{1,2},"TableName='[^']*'",'match'); 22 field = string(replace(field,["TableName='","'"," "],["","","_"])); 23 24 for k = 1:size(sbtab_excel{1,n},2) 25 26 if ~isempty(sbtab_excel{1,n}{2,k}) 27 subfield = sbtab_excel{1,n}{2,k}; 28 subfield =... 29 string(replace(subfield,["!",">",":"," "],["","","_","_"])); 30 31 sb.(field).(subfield)(:,1) = sbtab_excel{1,n}(3:end,k)'; 32 33 sb.(field).(subfield) = sb.(field).(subfield)... 34 (~cellfun('isempty', sb.(field).(subfield))); 35 end 36 end 37 end 38end 39end
Loads the SBtab saved in the .mat file and creates a MATLAB® struct that can be more easily parsed.
f_sbtab_to_model
Code
1function f_sbtab_to_model(stg,sb,mmf) 2% Saves the model in .mat, .sbproj and .xml format, while also creating a 3% file whith the data to run the model in all different experimental 4% settings defined in the sbtab 5 6modelobj = sbiomodel(stg.name); 7compObj = []; 8 9sbtab.species = cat(2,sb.Compound.Name,sb.Compound.InitialValue,... 10 sb.Compound.IsConstant,sb.Compound.Unit,sb.Compound.Location); 11 12sbtab.defpar = cat(2,sb.Parameter.Comment,sb.Parameter.Value_linspace,... 13 sb.Parameter.Unit); 14 15for n = 1:size(sb.Compartment.ID,2) 16 compObj{n} = addcompartment(modelobj, sb.Compartment.Name{n}); 17 set(compObj{n}, 'CapacityUnits', sb.Compartment.Unit{n}); 18 set(compObj{n}, 'Value', sb.Compartment.Size{n}); 19end 20 21for n = 1:size(sbtab.species,1) 22 23 for m = 1:size(compObj,2) 24 if string(compObj{m}.Name) == string(sb.Compound.Location{n}) 25 compartment_number_match = m; 26 end 27 end 28 29 addspecies (compObj{compartment_number_match}, sb.Compound.Name{n},sb.Compound.InitialValue{n}... 30 ,'InitialAmountUnits',sb.Compound.Unit{n}); 31end 32 33for n = 1:size(sbtab.defpar,1) 34 addparameter(modelobj,sb.Parameter.Name{n},... 35 sb.Parameter.Value_linspace{n},'ValueUnits',sb.Parameter.Unit{n},'Notes',sb.Parameter.Comment{n}); 36end 37 38for n = 1:size(sb.Reaction.ID,1) 39 40 if ischar(sb.Reaction.IsReversible{n}) 41 if contains(convertCharsToStrings(sb.Reaction.IsReversible{n}),"true") 42 reaction_name = strrep(sb.Reaction.ReactionFormula{n},'<=>',' <-> '); 43 else 44 reaction_name = strrep(sb.Reaction.ReactionFormula{n},'<=>',' -> '); 45 end 46 else 47 if sb.Reaction.IsReversible{n} 48 reaction_name = strrep(sb.Reaction.ReactionFormula{n},'<=>',' <-> '); 49 else 50 reaction_name = strrep(sb.Reaction.ReactionFormula{n},'<=>',' -> '); 51 end 52 end 53 54 reaction_name_compartment = reaction_name; 55 56 for m = 1:size(sb.Compound.Name,1) 57 reaction_name_compartment =... 58 insertBefore(string(reaction_name_compartment)," " +... 59 string(sb.Compound.Name{m})," " + string(sb.Reaction.Location{n})); 60 end 61 62 while contains(reaction_name_compartment,... 63 string(sb.Reaction.Location{n})+" "+string(sb.Reaction.Location{n})) 64 reaction_name_compartment =... 65 strrep(reaction_name_compartment,string(sb.Reaction.Location{n})+... 66 " "+string(sb.Reaction.Location{n})," "+sb.Reaction.Location{n}); 67 end 68 69 while contains(reaction_name_compartment," ") 70 reaction_name_compartment = strrep(reaction_name_compartment," "," "); 71 end 72 73 reaction_name_compartment = strrep(reaction_name_compartment,... 74 sb.Reaction.Location{n} + " ",sb.Reaction.Location{n}+"."); 75 reaction_name_compartment = string(sb.Reaction.Location{n})+"."+reaction_name_compartment; 76 77 reactionObj = addreaction(modelobj,reaction_name_compartment); 78 set(reactionObj,'ReactionRate',sb.Reaction.KineticLaw{n}); 79end 80 81for n = 1:size(sb.Compound.ID,1) 82 if ischar(sb.Compound.Assignment{n}) 83 if contains(convertCharsToStrings(sb.Compound.Assignment{n}),["true","True"]) 84 modelobj.species(n).BoundaryCondition = 1; 85 end 86 else 87 if sb.Compound.Assignment{n} == 1 88 modelobj.species(n).BoundaryCondition = 1; 89 end 90 end 91 if ischar(sb.Compound.Interpolation{n}) 92 if contains(convertCharsToStrings(sb.Compound.Interpolation{n}),["true","True"]) 93 modelobj.species(n).BoundaryCondition = 1; 94 end 95 else 96 if sb.Compound.Interpolation{n} == 1 97 modelobj.species(n).BoundaryCondition = 1; 98 end 99 end 100 if ischar(sb.Compound.IsConstant{n}) 101 if contains(convertCharsToStrings(sb.Compound.IsConstant{n}),["true","True"]) 102 modelobj.species(n).BoundaryCondition = 1; 103 end 104 else 105 if sb.Compound.IsConstant{n} == 1 106 modelobj.species(n).BoundaryCondition = 1; 107 end 108 end 109end 110 111sbtab.sim_time = [sb.Experiments.Sim_Time{:}]; 112 113species_INP_matcher ={}; 114for n = 1:size(sb.Compound.ID,1) 115 if isfield(sb.Experiments,"S"+(n-1)) 116 species_INP_matcher{size(species_INP_matcher,1)+1,1} = n; 117 end 118end 119 120for n = 1:size(sb.Experiments.ID,1) 121 startamount = cell(1,size(species_INP_matcher,1)); 122 nstartamount = 0; 123 nInputTime = 0; 124 nInput = 0; 125 nOutput = 0; 126 127 for m = 1:size(sb.Compound.ID,1) 128 if isfield(sb.Experiments,"S"+(m-1)) 129 nstartamount = nstartamount+1; 130 startamount{nstartamount} = eval("sb.Experiments.S"+(m-1)+"(n)"); 131 startAmountName(nstartamount) = sb.Compound.Name(m); 132 end 133 end 134 135 if isfield(sb.Experiments,"Normalize") 136 sbtab.datasets(n).Normalize = sb.Experiments.Normalize{n}; 137 else 138 sbtab.datasets(n).Normalize = []; 139 end 140 141 if isfield(eval(("sb.E")+(n-1)),"Time") 142 Data(n).Experiment.t = transpose(eval("[sb.E"+(n-1)+".Time{:}]")); 143 end 144 145 for m = 1:size(sb.Compound.ID,1) 146 if isfield(eval(("sb.E")+(n-1)+"I"),"Input_Time_S"+(m-1)) 147 nInputTime = nInputTime + 1; 148 sbtab.datasets(n).input_time{1,nInputTime} = ... 149 eval(("[sb.E") + (n-1) + "I.Input_Time_S" + (m-1) + "{:}]"); 150 end 151 if isfield(eval(("sb.E")+(n-1)+"I"),"S"+(m-1)) 152 nInput = nInput + 1; 153 sbtab.datasets(n).input_value{1,nInput} = ... 154 eval(("[sb.E") + (n-1) + "I.S" + (m-1) + "{:}]"); 155 sbtab.datasets(n).input{nInput} = char("S" + (m-1)); 156 end 157 end 158 159 for m = 1:size(sb.Output.ID,1) 160 if isfield(eval(("sb.E")+(n-1)),"Y"+(m-1)) 161 nOutput = nOutput+1; 162 Data(n).Experiment.x(:,nOutput) = eval(("[sb.E") +... 163 (n-1) + ".Y" + (m-1) + "{:}]"); 164 Data(n).Experiment.x_SD(:,nOutput) = eval(("[sb.E") +... 165 (n-1) + ".SD_Y" + (m-1) + "{:}]"); 166 sbtab.datasets(n).output{nOutput} = sb.Output.Name(m); 167% sbtab.datasets(n).output_value{nOutput} = ... 168% {convertStringsToChars(... 169% strrep(string(sb.Output.Location{m}) + "." +... 170% string(sb.Output.Name{m}) + " = " +... 171% string(sb.Output.Formula{m}),'eps','0.0001'))}; 172 sbtab.datasets(n).output_value{nOutput} = ... 173 {convertStringsToChars(... 174 string(sb.Output.Location{m}) + "." +... 175 string(sb.Output.Name{m}) + " = " +... 176 string(sb.Output.Formula{m}))}; 177 178 sbtab.datasets(n).output_unit{nOutput} = sb.Output.Unit{m}; 179 180 sbtab.datasets(n).output_name{nOutput} = ... 181 sb.Output.Name(m); 182 sbtab.datasets(n).output_ID{nOutput} = ... 183 sb.Output.ID(m); 184 sbtab.datasets(n).output_location{nOutput} = ... 185 sb.Output.Location(m); 186 end 187 end 188 sbtab.datasets(n).stg.outnumber = nOutput; 189 sbtab.datasets(n).start_amount = cat(2,startAmountName(:)... 190 ,transpose([startamount{:}]),species_INP_matcher); 191end 192 193if isfield(sb,"Expression") 194 for m = 1:size(sb.Expression.ID,1) 195 if isfield(sb.Expression,'Formula') 196 if isa(sb.Expression.Formula{m},'double') 197 addspecies (modelobj, char(sb.Expression.Name(m)),... 198 str2double(string(sb.Expression.Formula{m})),... 199 'InitialAmountUnits',sb.Expression.Unit{m}); 200 else 201 try 202 addspecies (modelobj, char(sb.Expression.Name(m)),0,... 203 'InitialAmountUnits',sb.Expression.Unit{m}); 204 catch 205 end 206 addrule(modelobj, char({convertStringsToChars(... 207 string(sb.Expression.Location{m}) + "." +... 208 string(sb.Expression.Name{m}) + " = " +... 209 string(sb.Expression.Formula{m}))}),... 210 'repeatedAssignment'); 211 end 212 else 213 addparameter(modelobj,char(sb.Expression.Name(m)),... 214 str2double(string(sb.Expression.DefaultValue{m})),... 215 'ValueUnits',sb.Expression.Unit{m}); 216 end 217 end 218end 219 220if isfield(sb,"Input") 221 for m = 1:size(sb.Input.ID,1) 222 if isfield(sb.Input,'Formula') 223 if isa(sb.Input.Formula{m},'double') 224 addspecies (modelobj, char(sb.Input.Name(m)),... 225 str2double(string(sb.Input.DefaultValue{m})),... 226 'InitialAmountUnits',sb.Input.Unit{m}); 227 else 228 try 229 addspecies (modelobj, char(sb.Input.Name(m)),0,... 230 'InitialAmountUnits',sb.Input.Unit{m}); 231 catch 232 end 233 addrule(modelobj, char({convertStringsToChars(... 234 string(sb.Input.Location{m}) + "." +... 235 string(sb.Input.Name{m}) + " = " +... 236 string(sb.Input.DefaultValue{m}))}),... 237 'repeatedAssignment'); 238 end 239 else 240 addparameter(modelobj,char(sb.Input.Name(m)),... 241 str2double(string(sb.Input.DefaultValue{m})),... 242 'ValueUnits',sb.Input.Unit{m}); 243 end 244 end 245end 246 247if isfield(sb,"Constant") 248 for m = 1:size(sb.Constant.ID,1) 249 addparameter(modelobj,char(sb.Constant.Name(m)),... 250 str2double(string(sb.Constant.Value{m})),... 251 'ValueUnits',sb.Constant.Unit{m}); 252 end 253end 254 255sbproj_model = mmf.model.data.sbproj_model; 256matlab_model = mmf.model.data.mat_model; 257data_model = mmf.model.data.data_model; 258xml_model = mmf.model.data.xml_model; 259 260sbiosaveproject(sbproj_model,'modelobj') 261 262save(matlab_model,'modelobj') 263 264save(data_model,'Data','sbtab','sb') 265 266sbmlexport(modelobj,xml_model) 267 268end
Reads information from the SBtab and saves the model in MATLAB (.mat, .sbproj) and SBML(.xml) format, while also creating a file whith the data to run the model in all different experimental settings defined in the SBtab.
Inputs - stg, sb
Saves - model file .mat, model file .sbproj, model file .xml, and data file
f_setup_input
Code
1function f_setup_input(stg,mmf) 2% Creates code that loads the inputs of each experiment into a .mat file, 3% and creates the code to read this inputs at runtime when the experiments 4% are being simulated, all this generated code is stored on the 5% Input_functions folder 6 7matlab_model = mmf.model.data.mat_model; 8data_model = mmf.model.data.data_model; 9inp_model_data = mmf.model.data.input_model_data; 10Model_folder = mmf.model.main; 11model_input = mmf.model.input_functions.input; 12 13%Find correct path for loading depending on the platform 14load(data_model,'sbtab') 15 16load(matlab_model,'modelobj'); 17 18for Exp_n = 1:size(sbtab.datasets,2) 19 20 for index = 1:size(sbtab.datasets(Exp_n).input,2) 21 if size(sbtab.datasets(Exp_n).input_value{index},2) > 100 22 23 input_name = strrep(modelobj.species(1+str2double(strrep(... 24 sbtab.datasets(Exp_n).input(index),'S',''))).name,".",""); 25 inpX = input_name + "X"; 26 inpT = input_name + "T"; 27 inph1 = input_name + "h1"; 28 inph2 = input_name + "h2"; 29 30 fullFileName = sprintf('%s.m',... 31 model_input + Exp_n + "_" + input_name ); 32 33 fileID = fopen(fullFileName, 'wt'); 34 35 inp_str = template1(); 36 inp_str = replace(inp_str,... 37 ["SBtab_name","Exp_n","input_name",... 38 "inpX", "inpT", "inph1", "inph2", "inp_model_data",... 39 "sbtab.sim_time(Exp_n)"],[stg.name,Exp_n,... 40 input_name, inpX, inpT, inph1,... 41 inph2, inp_model_data,... 42 sbtab.sim_time(Exp_n)]); 43 44 fprintf(fileID,inp_str); 45 fclose(fileID); 46 end 47 end 48end 49 50fullFileName = sprintf('%s.m',model_input + "_creator" ); 51 52fileID = fopen(fullFileName, 'wt'); 53fprintf(fileID, "function " + stg.name + "_input_creator(~)\n"); 54helper = 0; 55for Exp_n = 1:size(sbtab.datasets,2) 56 for index =1:size(sbtab.datasets(Exp_n).input,2) 57 if size(sbtab.datasets(Exp_n).input_value{index},2) > 100 58 input_name = strrep(modelobj.species(1+str2double(strrep(... 59 sbtab.datasets(Exp_n).input(index),'S',''))).name,".",""); 60 if helper == 0 61 helper = 1; 62 helper2 = Exp_n; 63 end 64 if index == 1 && Exp_n == helper2 65 fprintf(fileID,"load('" + data_model + "','sbtab');\n"); 66 inp_creator_str = template2(); 67 inp_creator_str = replace(inp_creator_str,... 68 ["Exp_n", "input_name", "index", "inp_model_data"],... 69 [Exp_n, input_name, index, inp_model_data]); 70 % fprintf(fileID,inp_creator_str); 71 else 72 inp_creator_str = template3(); 73 inp_creator_str = replace(inp_creator_str,... 74 ["Exp_n", "input_name", "index", "inp_model_data"],... 75 [Exp_n, input_name, index, inp_model_data]); 76 end 77 fprintf(fileID,inp_creator_str); 78 end 79 end 80end 81fprintf(fileID, "end\n"); 82fclose(fileID); 83 84addpath(genpath(Model_folder)); 85eval(stg.name + "_input_creator()"); 86end 87 88function inp_str = template1() 89inp_str =... 90 "function thisAmp = SBtab_name_inputExp_n_input_name(times)\n"+... 91 "persistent inpX\n"+... 92 "persistent inpT\n"+... 93 "persistent inph1\n"+... 94 "persistent inph2\n"+... 95 "if isempty(inpX)\n"+... 96 "Data = coder.load('inp_model_data','expExp_n_input_name');\n"+... 97 "inpX = Data.expExp_n_input_name(:,2);\n"+... 98 "inpT = Data.expExp_n_input_name(:,1);\n"+... 99 "inph1 = 1;\n"+... 100 "end\n"+... 101 "thisAmp = inpX(end);\n"+... 102 "if times ~= sbtab.sim_time(Exp_n)\n"+... 103 "if times == 0\n"+... 104 "inph1 = 1;\n"+... 105 "thisAmp = inpX(1);\n"+... 106 "else\n"+... 107 "while times > inpT(inph1)\n"+... 108 "inph1 = inph1 + 1;\n"+... 109 "end\n"+... 110 "while times < inpT(inph1-1)\n"+... 111 "inph1 = inph1-1;\n"+... 112 "end\n"+... 113 "inph2 = (inpT(inph1)-times)*1/(inpT(inph1)-inpT(inph1-1));\n"+... 114 "thisAmp = (inpX(inph1-1)*inph2 + inpX(inph1)*(1-inph2));\n"+... 115 "end\n"+... 116 "end\n"+... 117 "end"; 118end 119 120 121function inp_creator_str = template2() 122inp_creator_str = ... 123 "expExp_n_input_name(:,1) = sbtab.datasets(Exp_n).input_time{index};\n"+... 124 "expExp_n_input_name(:,2) = sbtab.datasets(Exp_n).input_value{index};\n"+... 125 "save('inp_model_data','expExp_n_input_name');\n"; 126end 127function inp_creator_str = template3() 128inp_creator_str = ... 129 "expExp_n_input_name(:,1) = sbtab.datasets(Exp_n).input_time{index};\n"+... 130 "expExp_n_input_name(:,2) = sbtab.datasets(Exp_n).input_value{index};\n"+... 131 "save('inp_model_data','expExp_n_input_name','-append');\n"; 132end
Creates code that loads the inputs of each experiment into a .mat file and the code to read these inputs at runtime when the experiments are being simulated. All this generated code is stored on the “Model/’model folder name’/Formulas” folder.
Inputs - stg
Saves - model-specific functions
f_build_model_exp
Code
1function f_build_model_exp(stg,sb,mmf) 2%Creates two .mat files for each experiment, with all the added rules, 3%species and parameters needed depending on the inputs and outputs 4%specified on the sbtab, one for the equilibrium simulation run and one for 5%the proper run 6 7data_model = mmf.model.data.data_model; 8mat_model = mmf.model.data.mat_model; 9model_exp_eq = mmf.model.data.model_exp.equilibration; 10model_exp_default = mmf.model.data.model_exp.default; 11model_exp_detail = mmf.model.data.model_exp.detail; 12 13%Find correct path for loading depending on the platform 14load(data_model,'Data','sbtab') 15load(mat_model,'modelobj'); 16 17model_run = cell(size(sb.Experiments.ID,1),1); 18configsetObj = cell(size(sb.Experiments.ID,1),1); 19 20for number_exp = 1:size(sb.Experiments.ID,1) 21 22 output_value = sbtab.datasets(number_exp).output_value; 23 output = sbtab.datasets(number_exp).output; 24 output_unit = sbtab.datasets(number_exp).output_unit; 25 input_time = sbtab.datasets(number_exp).input_time; 26 input_value = sbtab.datasets(number_exp).input_value; 27 input_species = sbtab.datasets(number_exp).input; 28 29 30 model_run{number_exp} = copyobj(modelobj); 31 configsetObj{number_exp} = getconfigset(model_run{number_exp}); 32 33 set(configsetObj{number_exp}, 'MaximumWallClock', stg.maxt); 34 set(configsetObj{number_exp}, 'StopTime', stg.eqt); 35 set(configsetObj{number_exp}.CompileOptions,... 36 'DimensionalAnalysis', stg.dimenanal); 37 set(configsetObj{number_exp}.CompileOptions,... 38 'UnitConversion', stg.UnitConversion); 39 set(configsetObj{number_exp}.SolverOptions,... 40 'AbsoluteToleranceScaling', stg.abstolscale); 41 set(configsetObj{number_exp}.SolverOptions,... 42 'RelativeTolerance', stg.reltol); 43 set(configsetObj{number_exp}.SolverOptions,... 44 'AbsoluteTolerance', stg.abstol); 45 set(configsetObj{number_exp}.SolverOptions, 'OutputTimes', stg.eqt); 46 set(configsetObj{number_exp}, 'TimeUnits', stg.simtime); 47 set(configsetObj{number_exp}.SolverOptions,... 48 'MaxStep', stg.maxstepeq); 49 50 set(configsetObj{number_exp}.SolverOptions,... 51 'AbsoluteToleranceStepSize', stg.abstolstepsize_eq); 52 53 54 model_exp = model_run{number_exp}; 55 config_exp = configsetObj{number_exp}; 56 57 save(model_exp_eq + number_exp + ".mat",'model_exp','config_exp') 58 59 sbiosaveproject(model_exp_eq + number_exp + ".sbproj",'model_exp') 60 61 set(configsetObj{number_exp}, 'StopTime', sbtab.sim_time(number_exp)); 62 63 set(configsetObj{number_exp}.SolverOptions, 'OutputTimes',... 64 Data(number_exp).Experiment.t); 65 66 set(configsetObj{number_exp}.SolverOptions, 'MaxStep', stg.maxstep); 67 68 for n = 1:size(output,2) 69 70 m = 0; 71 for k = 1:size(model_run{number_exp}.species,1) 72 if model_run{number_exp}.species(k).name == string(output{1,n}) 73 model_run{number_exp}.species(k).BoundaryCondition = 1; 74 m = 1; 75 end 76 end 77 78 if m == 0 79 if strcmp( output_unit{1,n},'dimensionless' ) 80 warning('off','SimBiology:InvalidSpeciesInitAmtUnits') 81 else 82 warning('on','SimBiology:InvalidSpeciesInitAmtUnits') 83 end 84 addspecies (model_run{number_exp}.Compartments(1),... 85 char(output{1,n}),0,... 86 'InitialAmountUnits',output_unit{1,n}); 87 end 88 89 addrule(model_run{number_exp}, char(output_value{1,n}),... 90 'repeatedAssignment'); 91 end 92 93 for j = 1:size(input_species,2) 94 95 input_indexcode = str2double(strrep(input_species(j),'S','')); 96 input_name = string(model_run{number_exp}.species(1 +... 97 input_indexcode).name); 98 99 if size(input_time{j},2) < 100 100 101 model_run{number_exp}.species(... 102 1 + input_indexcode).BoundaryCondition = 1; 103 for n = 1:size(input_time{j},2) 104 if ~isnan(input_time{j}(n)) 105 106 addparameter(model_run{number_exp},char("time_event_t_" + j + "_" + n),... 107 str2double(string(input_time{j}(n))),... 108 'ValueUnits',char(stg.simtime)); 109 110 addparameter(model_run{number_exp},char("time_event_r_" + j + "_" + n),... 111 str2double(string(input_value{j}(n))),... 112 'ValueUnits',char(model_run{number_exp}.species(1 +... 113 input_indexcode).InitialAmountUnits)); 114 115 addevent(model_run{number_exp}, ... 116 char("time>=time_event_t_" + j + "_" + n),... 117 cellstr(sbtab.datasets(number_exp).output_location{1} +... 118 "." + input_name + " = time_event_r_" + j + "_" + n)); 119 120 end 121 end 122 else 123 124 addrule(model_run{number_exp}, char(sbtab.datasets(... 125 number_exp).output_location{1} + "." + input_name +... 126 "=" + string(model_run{number_exp}.name)+ "_input" + ... 127 number_exp + "_" + input_name + "(time)"), 'repeatedAssignment'); 128 end 129 end 130 131 model_exp = model_run{number_exp}; 132 config_exp = configsetObj{number_exp}; 133 134 save(model_exp_default + number_exp + ".mat",'model_exp','config_exp') 135 136 sbiosaveproject(model_exp_default + number_exp + ".sbproj",'model_exp') 137 138 set(configsetObj{number_exp}.SolverOptions,'OutputTimes', []); 139 set(configsetObj{number_exp}.SolverOptions,'MaxStep', stg.maxstepdetail); 140 141 model_exp = model_run{number_exp}; 142 config_exp = configsetObj{number_exp}; 143 144 save(model_exp_detail + number_exp + ".mat",'model_exp','config_exp') 145 146end 147end
Creates two .mat files for each experiment, one for the equilibrium simulation run and one for the proper simulation. These files have all the added rules, species and parameters needed depending on the inputs and outputs specified on the SBtab.
Inputs - stg, sb
Saves - Ready to run models
Simulation and Scoring
f_sim_score
Code
1function [score,rst,rst_not_simd] = f_sim_score(parameters,stg,mmf) 2 3%Turn off Dimension analysis warning from simbiology 4warning('off','SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless') 5 6% Call the function that simulates the model 7rst = f_prep_sim(parameters,stg,mmf); 8 9% Call the function that scores 10rst = f_score(rst,stg,mmf); 11 12% Get the total score explicitly for optimization functions 13score = rst.st; 14 15rst_not_simd = rmfield( rst , 'simd'); 16end
Calls the function that runs the simulations and the function that scores the output of the runs.
Inputs
parameters - (double) Set of parameters that we are working on
Outputs
Calls - f_prep_sim, f_score
Loads
f_prep_sim
Code
1function rst = f_prep_sim(parameters,stg,mmf) 2 3% Save variables that need to be mantained over multiple function calls 4% persistent modelobj 5persistent sbtab 6persistent Data 7 8data_model = mmf.model.data.data_model; 9 10% Import the data on the first run 11if isempty(sbtab) 12 13 %Find correct path for loading depending on the platform 14 load(data_model,'Data','sbtab') 15end 16 17% Set the parameters that are going to be used for the simulation to the 18% default ones as definded in the SBTAB 19rt.par(:,1) = [sbtab.defpar{:,2}]; 20 21% Check if the parametrer needs to be set to the value relevant for Profile 22% Likelihood 23if isfield(stg,"PLind") && isfield(stg,"PLval") 24 parameters = [parameters(1:stg.PLind-1) stg.PLval parameters(stg.PLind:end)]; 25end 26 27% Iterate over all the parameters of the model 28for n = 1:size(rt.par,1) 29 30 % Check that a parameter should be changed from default 31 if stg.partest(n) > 0 32 33 % Set the parameters are being tested 34 rt.par(n) = 10.^(parameters(stg.partest(n,1))); 35 end 36 37 if isfield(stg,'tci') 38 39 % Check that there are thermodynamic constraints to implement 40 if ~isempty(stg.tci) 41 42 % Choose the parameters that need to be calculated with other 43 % parameters due to thermodynamic constraints 44 if ismember(n,stg.tci) 45 46 % Check that a parameter should be changed from default 47 if stg.partest(n) > 0 48 49 % Iterate over the parameters that need to be mutiplied 50 % for calculating the parameter that depends on the 51 % thermodynamic constraints 52 for m = 1:size(stg.tcm,2) 53 54 % Check that the parameter that is going to be used 55 % to calculate the parameter dependent on 56 % thermodynamic constraintsis is not the default 57 if stg.partest(stg.tcm(n,m),1) > 0 58 59 % Check if the parametrer needs to be set to 60 % the value relevant for Profile Likelihood 61 if isfield(stg,"PLind") 62 if stg.partest(stg.tcm(n,m),1) ==... 63 stg.PLind 64 parameters(stg.partest(... 65 stg.tcm(n,m),1))... 66 = stg.PLval; 67 end 68 end 69 70 % Make the appropriate multiplications to get 71 % the thermodinamicly constrained parameter 72 rt.par(n) = rt.par(n).*(10.^... 73 (parameters(stg.partest(... 74 stg.tcm(n,m),1)))); 75 else 76 77 % Make the appropriate multiplications to get 78 % the thermodinamicly constrained parameter 79 rt.par(n) = rt.par(n).*... 80 (sbtab.defpar{stg.tcm(n,m),2}); 81 end 82 end 83 84 % Iterate over the parameters that need to be divided 85 % for calculating the parameter that depends on the 86 % thermodynamic constraints 87 for m = 1:size(stg.tcd,2) 88 89 % Check that the parameter that is going to be used 90 % to calculate the parameter dependent on 91 % thermodynamic constraintsis is not the default 92 if stg.partest(stg.tcd(n,m),1) > 0 93 94 % Check if the parametrer needs to be set to 95 % the value relevant for Profile Likelihood 96 if isfield(stg,"PLind") 97 if stg.partest(stg.tcd(n,m),1) ==... 98 stg.PLind 99 parameters(stg.partest(... 100 stg.tcd(n,m),1))... 101 = stg.PLval; 102 end 103 end 104 % Make the appropriate divisions to get the 105 % thermodinamicly constrained parameter 106 rt.par(n) = rt.par(n)./(10.^... 107 (parameters(stg.partest(... 108 stg.tcd(n,m),1)))); 109 else 110 111 % Make the appropriate divisions to get the 112 % thermodinamicly constrained parameter 113 rt.par(n) = rt.par(n)./... 114 (sbtab.defpar{stg.tcd(n,m),2}); 115 end 116 end 117 end 118 end 119 end 120 end 121end 122 123% Set the start amount for the species in the model to 0 124rt.ssa = zeros(size(sbtab.species,1),max(stg.exprun)); 125 126% Initialize the results variable 127rst = []; 128rst.parameters = rt.par; 129% Iterate over all the experiments that are being run 130for n = stg.exprun 131 132 % Try catch used because iterations errors can happen unexectedly and 133 % we want to be able to continue simulations 134 try 135 136 % If the correct setting is chosen display messages to console 137 if stg.simcsl 138 disp("Running dataset number " + n +... 139 " of " + stg.exprun(end)) 140 end 141 142 % Check that this is not the first time the experiments are being 143 % run and that the start values for the species are different from 144 % the previous experiment 145 if n ~= stg.exprun(1) && ... 146 min([sbtab.datasets(n).start_amount{:,2}] ==... 147 [sbtab.datasets(max(1,stg.exprun(... 148 find(stg.exprun==n)-1))).start_amount{:,2}]) 149 150 151 % Set the values of the start amounts to the values obtained 152 % after the first equilibration 153 rt.ssa(:,n) =... 154 rt.ssa(:,stg.exprun(find(stg.exprun==n)-1)); 155 if stg.simdetail 156 rt.ssa(:,n+2*stg.expn) =... 157 rt.ssa(:,stg.exprun(find(stg.exprun==n)-1)); 158 end 159 else 160 161 % Iterate over the numbre of species that need a starting value 162 % different than 0 163 for j = 1:size(sbtab.datasets(n).start_amount,1) 164 165 % Set the start amount of the species to the number defined 166 % in the sbtab for each experiment 167 rt.ssa(sbtab.datasets(n... 168 ).start_amount{j,3},n+stg.expn) =... 169 sbtab.datasets(n).start_amount{j,2}; 170 end 171 172 % Equilibrate the model 173 rst = f_sim(n+stg.expn,stg,rt,rst,mmf); 174 175 for j = 1:size(sbtab.species,1) 176 177 % Set the starting amount for species that after 178 % equilibrium have very low values to zero 179 if rst.simd{n+stg.expn}.Data(end,j) < 1.0e-15 180 rt.ssa(j,n) = 0; 181 182 % Set the starting amount for the rest of the species 183 else 184 rt.ssa(j,n) =... 185 rst.simd{n+stg.expn}.Data(end,j); 186 if stg.simdetail 187 rt.ssa(j,n+2*stg.expn) =... 188 rst.simd{n+stg.expn}.Data(end,j); 189 end 190 end 191 end 192 end 193 194 % Simulate the model 195 rst = f_sim(n,stg,rt,rst,mmf); 196 197 try 198 if stg.simdetail 199 rst = f_sim(n+2*stg.expn,stg,rt,rst,mmf); 200 end 201 catch 202 end 203 204 % Check If the times of the simultaion output and the simulation 205 % data from SBTAB match, if they don't it means that the simulator 206 % didn't had enough time to run the model (happens in some 207 % unfavorable configurations of parameters, controlled by stg.maxt 208 if size(Data(n).Experiment.t,1) ~=... 209 size(rst.simd{n}.Data(:,end),1) 210 211 % Set the simulation output to be 0, this is a non function 212 % value that the score function expects in simulations that did 213 % not worked properly 214 rst.simd{n} = 0; 215 end 216 catch 217 218 % Set the simulation output to be 0, this is a non function value 219 % that the score function expects in simulations that did not 220 % worked properly 221 rst.simd{n} = 0; 222 end 223end 224end
Prepares the simulation making sure that an equilibration is preformed when necessary before running the main simulation.
Inputs
stg - stg.name, stg.partest, stg.tci, stg.tcm, stg.tcd, stg.exprun, stg.simcsl, stg.expn
parameters - (double) Set of parameters that we are working on
Created Variables
rt
rt.ssa - (double) steady state amounts
rt.par - (double) All parameters of the model, takes the default ones from SBtab and then replaces the ones being worked on.
Outputs
Calls - f_sim
f_sim
Code
1function rst = f_sim(exp_n,stg,rt,rst,mmf) 2 3% Save variables that need to be mantained over multiple function calls 4persistent model_run 5persistent config_run 6 7% If the function is called for the first time, load the appropriate model 8% and compile the code for simulation run 9if isempty(model_run) 10 11 warning('off','SimBiology:InvalidSpeciesInitAmtUnits') 12 13 %Generate an empty array to be populated with the model suited for each 14 %equilibration and experiment% 15 model_run = cell(1,stg.expn*2); 16 config_run = cell(1,stg.expn*2); 17 18 model_exp_default = mmf.model.data.model_exp.default; 19 model_exp_eq = mmf.model.data.model_exp.equilibration; 20 model_exp_detail = mmf.model.data.model_exp.detail; 21 22 % Iterate over the experiments thar are being run 23 for n = stg.exprun 24 25 if stg.simdetail 26 % Load the models for equilibrium 27 load(model_exp_detail + n + ".mat",'model_exp','config_exp') 28 29 % Place the loaded models in the correct place in the array, 30 % models for equilibrium are set to be on the second half of 31 % the array 32 model_run{n+2*stg.expn} = model_exp; 33 config_run{n+2*stg.expn} = config_exp; 34 35 % Compile the matlab code that is going to simulate the model 36 % using matlab built in function if the option is selected in 37 % settings 38 if stg.sbioacc 39 sbioaccelerate(model_run{n+2*stg.expn},config_run{n+2*stg.expn}); 40 end 41 end 42 % Load the models for equilibrium 43 load(model_exp_eq + n + ".mat",'model_exp','config_exp') 44 45 % Place the loaded models in the correct place in the array, models 46 % for equilibrium are set to be on the second half of the array 47 model_run{n+stg.expn} = model_exp; 48 config_run{n+stg.expn} = config_exp; 49 50 % Compile the matlab code that is going to simulate the model using 51 % matlab built in function if the option is selected in settings 52 if stg.sbioacc 53 sbioaccelerate(model_run{n+stg.expn},config_run{n+stg.expn}); 54 end 55 56 % Load the models for main run 57 load(model_exp_default + n + ".mat",'model_exp','config_exp') 58 59 % Place the loaded models in the correct place in the array, models 60 % for main run are set to be on the first half of the array 61 model_run{n} = model_exp; 62 config_run{n} = config_exp; 63 64 % Compile the matlab code that is going to simulate the model using 65 % matlab built in function if the option is selected in settings 66 if stg.sbioacc 67 sbioaccelerate(model_run{n},config_run{n}); 68 end 69 end 70end 71 72% substitute the start amount of the species in the model with the correct 73% ones for simulations 74set(model_run{exp_n}.species(1:size(rt.ssa(:,exp_n),1)),{'InitialAmount'},num2cell(rt.ssa(:,exp_n))) 75 76% Substitute the values of the parameters in the model for the correct one 77% for simultaions 78set(model_run{exp_n}.parameters(1:size(rt.par,1)),{'Value'},num2cell(rt.par)) 79 80%simulate the model using matlab built in function 81rst.simd{exp_n} = sbiosimulate(model_run{exp_n},config_run{exp_n}); 82end
Simulates the model with the provided configurations. The first time it is run it loads a representation of the model and the simulation, and compiles this information to C code.
Inputs
exp_n - (double) Unique number to identify the model for each experiment or equilibrium reaction (it needs a new model object for each one)
rt
rt.ssa - (double) steady state amounts
rt.par - (double) All parameters of the model, takes the default ones from SBtab and then replaces the ones being worked on.
Outputs
Calls - Sbioaccelerate, Sbiosimulate
Loads - Ready to run model, Ready to run model equilibration
f_score
Code
1function rst = f_score(rst,stg,mmf) 2 3persistent sbtab 4persistent Data 5 6data_model = mmf.model.data.data_model; 7 8% Import the data on the first run 9if isempty(Data) 10 load(data_model,'Data','sbtab') 11end 12 13% Iterate over the number of experiments 14for n = stg.exprun 15 16 % Iterate over the number of datasets per experiment 17 for j = 1:size(sbtab.datasets(n).output,2) 18 19 % Calculate score per dataset if there are no errors 20 if rst.simd{n} ~= 0 21 22 data = Data(n).Experiment.x(:,j); 23 data_sd = Data(n).Experiment.x_SD(:,j); 24 number_points = size(Data(n).Experiment.x(:,j),1); 25 sim_results = f_normalize(rst,stg,n,j,mmf); 26 rst.xfinal{n,1}(j) = sim_results(end); 27 28 % Calculate score using formula that accounts for normalization 29 % with the starting point of the result 30 if stg.useLog == 4 31 rst.sd(j,n) = sum(((data-sim_results)./... 32 (data_sd*sqrt(number_points))).^2); 33 rst.sdtest(j,n) = sum(((data-sim_results)./... 34 (data_sd*sqrt(number_points))).^2); 35 36% elseif stg.useLog == 5 37% rst.sd{n,1}(j) = sum(((data-sim_results)./... 38% (data_sd)).^2); 39 else 40 rst.sd(j,n) = sum(((data-sim_results)./... 41 (data_sd)).^2)/(number_points); 42 rst.sdtest(j,n) = sum(((data-sim_results)./... 43 (data_sd)).^2)/(number_points); 44 end 45 46 % If there are errors output a very high score value (10^10) 47 elseif rst.simd{n} == 0 || rst.sd(n,j) == inf 48 49 rst.sd(n,j) = stg.errorscore; 50 rst.xfinal{n,1}(j) = 0; 51% rst.sdtest(j,n) = stg.errorscore; 52 end 53% rst.sd{n,1}(j) 54% rst.sdtest(j,n) 55% max(0,log10(rst.sd{n,1}(j))); 56% max(0,log10(rst.sdtest(j,n))); 57 58 % Calculate the log10 of dataset score if option selected 59 if stg.useLog == 1 60 rst.sd(j,n) = max(0,log10(rst.sd(j,n))); 61 end 62 end 63 64 % Calculate score per experiment 65 rst.se(n,1) = sum(rst.sd(:,n)); 66% rst.setest(n,1) = sum(rst.sdtest(:,n)); 67 68 % Calculate the log10 of experiment score if option selected 69 if stg.useLog == 2 70 rst.se(n,1) = log10(rst.se(n,1)); 71 end 72end 73 74% Calculate score per experiment 75rst.st = sum(rst.se); 76 77% Calculate the log10 of total score if option selected 78if stg.useLog == 3 79 rst.st = log10(rst.st); 80end 81end
Uses the results from the simulation of the model and the Data provided via the SBTAB to calculate a score for a given parameter set.
Analysis
f_analysis
Code
1function rst = f_analysis(stg,analysis,mmf,analysis_options) 2if contains(analysis,analysis_options(1)) 3 disp("Starting " + analysis_options(1)) 4 rst.diag = f_diagnostics(stg,mmf); 5 disp(analysis_options(1) + " completed successfully") 6end 7 8if contains(analysis,analysis_options(2)) 9 disp("Starting " + analysis_options(2)) 10 rst.opt = f_opt(stg,mmf); 11 disp(analysis_options(2) + " completed successfully") 12end 13 14if contains(analysis,analysis_options(3)) 15 disp("Starting " + analysis_options(3)) 16 rst.gsa = f_gsa(stg,mmf); 17 disp(analysis_options(3) + " completed successfully") 18end 19 20if contains(analysis,analysis_options(4)) 21 disp("Starting " + analysis_options(4)) 22 rst.PLA = f_PL_m(stg,mmf); 23 disp(analysis_options(4) + " completed successfully") 24end 25end
Calls the proper analysis functions depending on the analysis that was chosen on the settings file. The supported analysis right now are:
Inputs
analysis - (string) analysis being run (stg.analysis)
Outputs - rst
Diagnostics
Code
1function rst = f_diagnostics(stg,mmf) 2 3% Run the model and obtain scores for fitness Multi Core 4if stg.optmc 5 disp("Running the model and obtaining Scores (Multicore)") 6 7 pa = stg.pa; 8 % Iterate over the parameter arrays to be tested 9 parfor n = stg.pat 10 [~,rst(n),~] = f_sim_score(pa(n,:),stg,mmf); 11 end 12 13 % Run the model and obtain scores for fitness single Core 14else 15 disp("Running the model and obtaining Scores (Singlecore)") 16 17 % Iterate over the parameter arrays to be tested 18 for n = stg.pat 19 [~,rst(n),~] = f_sim_score(stg.pa(n,:),stg,mmf); 20 end 21end 22end
Inputs
Outputs - rst (diagnostics results)
Optimization
Code
1function rst = f_opt(stg,mmf) 2% Call function to run fmincon optimization algorithm if chosen in settings 3 4if stg.fmincon 5 rst(1) = f_opt_fmincon(stg,mmf); 6end 7 8% Call function to run simulated annealing optimization algorithm if chosen 9% in settings 10if stg.sa 11 rst(2) = f_opt_sa(stg,mmf); 12end 13 14% Call function to run pattern search optimization algorithm if chosen in 15% settings 16if stg.psearch 17 rst(3) = f_opt_psearch(stg,mmf); 18end 19 20% Call function to run genetic algorithm optimization if chosen in settings 21if stg.ga 22 rst(4) = f_opt_ga(stg,mmf); 23end 24 25% Call function to run Particle swarm optimization algorithm if chosen in 26% settings 27if stg.pswarm 28 rst(5) = f_opt_pswarm(stg,mmf); 29end 30 31% Call function to run Surrogate optimization algorithm if chosen in 32% settings 33if stg.sopt 34 rst(6) = f_opt_sopt(stg,mmf); 35end 36end
Calls the correct optmizer or optimizers that have been chosen in the settings file.
Inputs
Outputs - rst (optimization results)
Code
1function [spoint,spop] = f_opt_start(stg) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Optimization Start method 1 7if stg.osm == 1 8 9 % Get a random starting point or group of starting points, if using 10 % multistart, inside the bounds 11 spoint = lhsdesign(stg.msts,stg.parnum).*(stg.ub-stg.lb)+stg.lb; 12 13 % Get a group of ramdom starting points inside the bounds 14 spop = lhsdesign(stg.popsize,stg.parnum).*(stg.ub-stg.lb)+stg.lb; 15 16 % Optimization Start method 2 17elseif stg.osm == 2 18 19 % Get a random starting point or group of starting points, if using 20 % multistart, near the best point 21 spoint = stg.bestpa - stg.dbpa +... 22 (stg.dbpa*2*lhsdesign(stg.msts,stg.parnum)); 23 24 % Get a group of ramdom starting points near the best point 25 spop = stg.bestpa - stg.dbpa +... 26 (stg.dbpa*2*lhsdesign(stg.popsize,stg.parnum)); 27end 28end
Inputs
Outputs
spoint - (double) starting parameter set for the optimization
spop - (double) Starting parameter sets for multiple start optimizations
Code
f_opt_fmincon
1function rst = f_opt_fmincon(stg,mmf) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Set the starting point for the optimization 7[startpoint,~] = f_opt_start(stg); 8 9% Get the optimization options from settings 10options = stg.fm_options; 11options.UseParallel = stg.optmc; 12 13% Display console messages if chosen in settings 14if stg.optcsl 15 options.Display = 'iter-detailed'; 16end 17 18% Display plots if chosen in settings 19if stg.optplots 20 options.PlotFcn = ... 21 {@optimplotx,@optimplotfunccount,@optimplotfval,... 22 @optimplotstepsize,@optimplotfirstorderopt}; 23end 24 25% Optimize the model with multiple starting points if chosen in settings 26if stg.mst 27 parfor n = 1:stg.msts 28 disp(string(n)) 29 [x(n,:),fval(n),exitflag(n),output(n)] =... 30 fmincon(@(x)f_sim_score(x,stg,mmf),startpoint(n,:),... 31 [],[],[],[],stg.lb,stg.ub,[],options); 32 end 33 34 % Optimize the model 35else 36 [x(1,:),fval(1),exitflag(1),output(1)] =... 37 fmincon(@(x)f_sim_score(x,stg,mmf),(stg.lb+stg.ub)/2,... 38 [],[],[],[],stg.lb,stg.ub,[],options); 39end 40 41% Save results 42rst.name = 'fmincon'; 43rst.x = x; 44rst.fval = fval; 45rst.exitflag = exitflag; 46rst.output = output; 47endf_opt_sa
1function rst = f_opt_sa(stg,mmf) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Set the starting point for the optimization 7[startpoint,~] = f_opt_start(stg); 8 9% Get the optimization options from settings 10options = stg.sa_options; 11options.MaxTime = stg.optt; 12options.InitialTemperature = ones(1,stg.parnum)*2; 13 14% Display console messages if chosen in settings 15if stg.optcsl 16 options.Display = 'iter'; 17end 18 19% Display plots if chosen in settings 20if stg.optplots 21 options.PlotFcn =... 22 {@saplotbestf,@saplottemperature,@saplotf,@saplotstopping,@saplotx}; 23end 24 25% Optimize the model with multiple starting points if chosen in settings 26if stg.mst 27 parfor n = 1:stg.msts 28 disp(string(n)) 29 [x(n,:),fval(n),exitflag(n),output(n)] =... 30 simulannealbnd(@(x)f_sim_score(x,stg,mmf),startpoint(n,:),... 31 stg.lb,stg.ub,options); 32 end 33 34 % Optimize the model 35else 36 [x(1,:),fval(1),exitflag(1),output(1)] =... 37 simulannealbnd(@(x)f_sim_score(x,stg,mmf),(stg.lb+stg.ub)/2,... 38 stg.lb,stg.ub,options); 39end 40 41% Save results 42rst.name = 'Simulated annealing'; 43rst.x = x; 44rst.fval = fval; 45rst.exitflag = exitflag; 46rst.output = output; 47endf_opt_psearch
1function rst = f_opt_psearch(stg,mmf) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Set the starting point for the optimization 7[startpoint,~] = f_opt_start(stg); 8 9% Get the optimization options from settings 10options = stg.psearch_options; 11options.MaxTime = stg.optt; 12options.UseParallel = stg.optmc; 13 14% Display console messages if chosen in settings 15if stg.optcsl 16 options.Display = 'iter'; 17end 18 19% Display plots if chosen in settings 20if stg.optplots 21 options.PlotFcn = ... 22 {@psplotbestf,@psplotfuncount,@psplotmeshsize,@psplotbestx}; 23end 24 25% Optimize the model with multiple starting points if chosen in settings 26if stg.mst 27 parfor n = 1:stg.msts 28 disp(string(n)) 29 [x(n,:),fval(n),exitflag(n),output(n)] =... 30 patternsearch(@(x)f_sim_score(x,stg,mmf),startpoint(1,:),[],[], ... 31 [],[],stg.lb,stg.ub,[],options); 32 end 33 34 % Optimize the model 35else 36 [x(1,:),fval(1),exitflag(1),output(1)] =... 37 patternsearch(@(x)f_sim_score(x,stg,mmf),(stg.lb+stg.ub)/2,[],[], ... 38 [],[],stg.lb,stg.ub,[],options); 39 40end 41 42% Save results 43rst.name = 'Pattern search'; 44rst.x = x; 45rst.fval = fval; 46rst.exitflag = exitflag; 47rst.output = output; 48endf_opt_ga
1function rst = f_opt_ga(stg,mmf) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Get the optimization options from settings 7options = stg.ga_options; 8options.MaxTime = stg.optt; 9options.UseParallel = stg.optmc; 10options.PopulationSize = stg.popsize; 11 12% Display console messages if chosen in settings 13if stg.optcsl 14 options.Display = 'iter'; 15end 16 17% Display plots if chosen in settings 18if stg.optplots 19 options.PlotFcn = {@gaplotbestf,... 20 @gaplotexpectation,@gaplotrange,@gaplotscorediversity,... 21 @gaplotstopping,@gaplotscores,@gaplotdistance,... 22 @gaplotselection}; 23end 24 25% Set the starting population for the optimization 26[~,startpop] = f_opt_start(stg); 27options.InitialPopulationMatrix = startpop; 28 29% Optimize the model with multiple starting points if chosen in settings 30if stg.mst 31 parfor n = 1:stg.msts 32 disp(string(n)) 33 [x(n,:),fval(n)] = ga(@(x)f_sim_score(x,stg,mmf),stg.parnum,... 34 [],[],[],[],stg.lb,stg.ub,[],options); 35 end 36 37 % Optimize the model 38else 39 [x(1,:),fval(1)] = ga(@(x)f_sim_score(x,stg,mmf),stg.parnum,... 40 [],[],[],[],stg.lb,stg.ub,[],options); 41end 42 43% Save results 44rst.name = 'Genetic algorithm'; 45rst.x = x; 46rst.fval = fval; 47rst.exitflag = []; 48rst.output = []; 49endf_opt_pswarm
1function rst = f_opt_pswarm(stg,mmf) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Get the optimization options from settings 7options = stg.pswarm_options; 8options.MaxTime = stg.optt; 9options.UseParallel = stg.optmc; 10options.SwarmSize = stg.popsize; 11 12% Display console messages if chosen in settings 13if stg.optcsl 14 options.Display = 'iter'; 15end 16 17% Display plots if chosen in settings 18if stg.optplots 19 options.PlotFcn = @pswplotbestf; 20end 21 22% Set the starting population for the optimization 23[~,startpop] = f_opt_start(stg); 24options.InitialSwarmMatrix = startpop; 25 26% Optimize the model with multiple starting points if chosen in settings 27if stg.mst 28 parfor n = 1:stg.msts 29 disp(string(n)) 30 [x(n,:),fval(n),exitflag(n),output(n)] =... 31 particleswarm(@(x)f_sim_score(x,stg,mmf),... 32 stg.parnum,stg.lb,stg.ub,options); 33 end 34 35 % Optimize the model 36else 37 [x(1,:),fval(1),exitflag(1),output(1)] =... 38 particleswarm(@(x)f_sim_score(x,stg,mmf),... 39 stg.parnum,stg.lb,stg.ub,options); 40end 41 42% Save results 43rst.name = 'Particle swarm'; 44rst.x = x; 45rst.fval = fval; 46rst.exitflag = exitflag; 47rst.output = output; 48endf_opt_sopt
1function rst = f_opt_sopt(stg,mmf) 2 3% Set the randomm seed for reproducibility 4rng(stg.rseed); 5 6% Get the optimization options from settings 7options = stg.sopt_options; 8options.MaxTime = stg.optt; 9options.UseParallel = stg.optmc; 10 11% Display console messages if chosen in settings 12if stg.optcsl 13 options.Display = 'iter'; 14end 15 16% Display plots if chosen in settings 17if stg.optplots 18 options.PlotFcn = @surrogateoptplot; 19end 20 21% Set the starting population for the optimization 22[~,startpop] = f_opt_start(stg); 23options.InitialPoints = startpop; 24 25% Optimize the model with multiple starting points if chosen in settings 26if stg.mst 27 parfor n = 1:stg.msts 28 disp(string(n)) 29 [x(n,:),fval(n),exitflag(n),output(n)] =... 30 surrogateopt(@(x)f_sim_score(x,stg,mmf),stg.lb,stg.ub,options); 31 end 32 33 % Optimize the model 34else 35 [x(1,:),fval(1),exitflag(1),output(1)] =... 36 surrogateopt(@(x)f_sim_score(x,stg,mmf),stg.lb,stg.ub,options); 37end 38 39% Save results 40rst.name = 'Surrogate optimization'; 41rst.x = x; 42rst.fval = fval; 43rst.exitflag = exitflag; 44rst.output = output; 45end
These functions call built in MATLAB® functions that perform parameter optimization . For furher information relating to how these optimizers work please follow the links to the MATLAB® documentation. Optimizers used:
f_opt_fmincon - fmincon
f_opt_sa - Simmulated annealing
f_opt_psearch - Pattern search
f_opt_ga - Genetic algorihtm
f_opt_pswarm - Particle swarm
f_opt_sopt - Surrogate optmization
Inputs - stg
Outputs - Optimization results
Global Sensitivity Analysis
Code
1function rst = f_gsa(stg,mmf) 2 3rst = f_make_par_samples(stg); 4 5rst = f_make_output_sample(rst,stg,mmf); 6 7rst = f_calc_sensitivities(rst,stg); 8end
Calls the global sensitivity analysis functions in the correct order.
Code
1function rst= f_make_par_samples(stg) 2% Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J. 3% comp. neuroscience 27.3 (2009): 471.) 4 5% MAKE SAMPLE MATRICES 6M1 = zeros(stg.sansamples, stg.parnum); % Pre-allocate memory for data 7M2 = zeros(stg.sansamples, stg.parnum); 8N = zeros(stg.sansamples, stg.parnum, stg.parnum); 9rng(stg.rseed) 10 11% Create a distribution for each parameter acording to settings(Note that 12% the sampling is being performed in log space as the parameters and its 13% bounds are in log space) 14for i=1:stg.parnum 15 % Uniform distribution truncated at the parameter bounds 16 if stg.sasamplemode == 0 17 M1(:,i) = stg.lb(i) +... 18 (stg.ub(i)-stg.lb(i)).*rand(1,stg.sansamples); 19 M2(:,i) = stg.lb(i) +... 20 (stg.ub(i)-stg.lb(i)).*rand(1,stg.sansamples); 21 % Normal distribution with mu as the best value for a parameter and 22 % sigma as stg.sasamplesigma truncated at the parameter bounds 23 elseif stg.sasamplemode == 1 24 pd(i) = makedist('Normal','mu',stg.bestpa(i),... 25 'sigma',stg.sasamplesigma); 26 t(i) = truncate(pd(i),stg.lb(i),stg.ub(i)); 27 r{i} = random(t(i),stg.sansamples,1); 28 r2{i} = random(t(i),stg.sansamples,1); 29 M1(:,i) = r{i}; 30 M2(:,i) = r2{i}; 31 % Same as 1 without truncation 32 elseif stg.sasamplemode == 2 33 pd(i) = makedist('Normal','mu',stg.bestpa(i),... 34 'sigma',stg.sasamplesigma); 35 r{i} = random(pd(i),stg.sansamples,1); 36 r2{i} = random(pd(i),stg.sansamples,1); 37 M1(:,i) = r{i}; 38 M2(:,i) = r2{i}; 39 % Normal distribution centered at the mean of the parameter bounds 40 % and sigma as stg.sasamplesigma truncated at the parameter bounds 41 elseif stg.sasamplemode == 3 42 pd(i) = makedist('Normal','mu',... 43 stg.lb(i) + (stg.ub(i)-stg.lb(i))/2,'sigma',stg.sasamplesigma); 44 t(i) = truncate(pd(i),stg.lb(i),stg.ub(i)); 45 r{i} = random(t(i),stg.sansamples,1); 46 r2{i} = random(t(i),stg.sansamples,1); 47 M1(:,i) = r{i}; 48 M2(:,i) = r2{i}; 49 % Same as 3 without truncation. 50 elseif stg.sasamplemode == 4 51 pd(i) = makedist('Normal','mu',... 52 stg.lb(i) + (stg.ub(i)-stg.lb(i))/2,'sigma',stg.sasamplesigma); 53 r{i} = random(pd(i),stg.sansamples,1); 54 r2{i} = random(pd(i),stg.sansamples,1); 55 M1(:,i) = r{i}; 56 M2(:,i) = r2{i}; 57 end 58end 59 60for i=1:stg.parnum 61 % Replace the i:th column in M2 by the i:th column from M1 to obtain Ni 62 N(:,:,i) = M2; 63 N(:,i,i) = M1(:,i); 64end 65 66rst.M1=M1; 67rst.M2=M2; 68rst.N=N;
Creates parameter sets samples with specific parameter distributions that are used to perform the global sensitivity analysis.
Inputs
stg - stg.sansamples, stg.parnum, stg.sasamplemode, stg.ub, stg.lb
Code inspired by Geir Halnes et al. 2009 paper.
Code
1function rst = f_make_output_sample(rst,stg,mmf) 2% Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J. 3% comp. neuroscience 27.3 (2009): 471.) 4 5nSamples = stg.sansamples; 6nPars = stg.parnum; 7parameter_array = zeros(nSamples,nPars); 8progress = 1; 9time_begin = datetime; 10D = parallel.pool.DataQueue; 11 12afterEach(D, @progress_track); 13 14for i=1:nSamples 15 parameter_array(i,:) = rst.M1(i,:); 16end 17 18parfor i=1:nSamples 19 [~,~,RM1(i)] = f_sim_score(parameter_array(i,:),stg,mmf); 20 send(D, "GSA M1 "); 21end 22disp("GSA M1 Runtime: " + string(datetime - time_begin) +... 23 " All " + nSamples + " samples executed") 24 25% [RM1(i).sd{:}] 26% RM1(i).sdtest(:,:) 27% reshape(RM1(i).sd(:,:),1,[]) 28 29for i=1:nSamples 30 fM1.sd(i,:) = reshape(RM1(i).sd(:,:),1,[]); 31 fM1.se(i,:) = RM1(i).se(:); 32 fM1.st(i,:) = RM1(i).st; 33 fM1.xfinal(i,:) = [RM1(i).xfinal{:}]; 34end 35 36rst.fM1 = fM1; 37clear a FM1 38 39for i=1:nSamples 40 parameter_array(i,:)= rst.M2(i,:); 41end 42 43progress = 1; 44parfor i=1:nSamples 45 [~,~,RM2(i)] = f_sim_score(parameter_array(i,:),stg,mmf); 46 send(D, "GSA M2 "); 47end 48disp("GSA M2 Runtime: " + string(datetime - time_begin) +... 49 " All " + nSamples + " samples executed") 50 51for i=1:nSamples 52 fM2.sd(i,:) = reshape(RM2(i).sd(:,:),1,[]); 53 fM2.se(i,:) = RM2(i).se(:); 54 fM2.st(i,:) = RM2(i).st; 55 fM2.xfinal(i,:) = [RM2(i).xfinal{:}]; 56end 57 58rst.fM2 = fM2; 59clear b FM2 60 61for i=1:nSamples 62 for j=1:nPars 63 parameter_array(i,:,j)= rst.N(i,:,j); 64 end 65end 66 67progress = 1; 68parfor i=1:nSamples 69 for j=1:nPars 70 [~,~,RN{i,j}] = f_sim_score(parameter_array(i,:,j),stg,mmf); 71 end 72 send(D, "GSA N "); 73end 74disp("GSA N Runtime: " + string(datetime - time_begin) +... 75 " All " + nSamples + " samples executed") 76 77for i=1:nSamples 78 for j=1:nPars 79 fN.sd(i,:,j) = reshape(RN{i,j}.sd(:,:),1,[]); 80 fN.se(i,:,j) = RN{i,j}.se(:); 81 fN.st(i,:,j) = RN{i,j}.st; 82 fN.xfinal(i,:,j) = [RN{i,j}.xfinal{:}]; 83 end 84end 85 86rst.fN = fN; 87clear c FN 88 89 function progress_track(name) 90 progress = progress + 1; 91 if mod(progress,ceil(nSamples/10)) == 0 && progress ~= nSamples 92 disp(name + "Runtime: " + string(datetime - time_begin) +... 93 " Samples:" + progress + "/" + nSamples) 94 end 95 end 96end
For each parameter set given in the matrices M1, M2, and N it runs the function f_sim_score generating new matrices fM1, fM2, and fN respectively.
Inputs - M1, M2, N, stg.sansamples, stg.parnum,
Code inspired by Geir Halnes et al. 2009 paper.
Code
1function rst = f_calc_sensitivities(rst,stg) 2 3rst = remove_sim_error(rst,stg); 4 5[rst.SiQ,rst.SiTQ,rst.Si,rst.SiT] = bootstrap(rst,stg); 6end 7 8function [SiQ,SiTQ,Si,SiT]=bootstrap(rst,stg) 9% calculates confidence intervals. 10fM1 = rst.fM1; 11fM2 = rst.fM2; 12fN = rst.fN; 13 14scores_names_list = ["sd","se","st","xfinal"]; 15 16[Si,SiT] = bootstrap_h(fM1,fM2,fN,stg,scores_names_list); 17 18if (isempty(stg.gsabootstrapsize)) 19 stg.gsabootstrapsize=ceil(sqrt(size(fM1.sd))); 20end%if 21 22fM1q = cell(stg.gsabootstrapsize,1); 23fM2q = cell(stg.gsabootstrapsize,1); 24fNq = cell(stg.gsabootstrapsize,1); 25 26parfor j=1:stg.gsabootstrapsize 27 [SiQh{j},SiTQh{j}] = bootstrap_hq(fM1,fM2,fN,stg,scores_names_list,j); 28end%parfor 29 30for n = 1:size(scores_names_list,2) 31 for j=1:stg.gsabootstrapsize 32 eval("SiQ." + scores_names_list(n) + "(j,:,:) = SiQh{j}." + scores_names_list(n) + ";") 33 eval("SiTQ." + scores_names_list(n) + "(j,:,:) = SiTQh{j}." + scores_names_list(n) + ";") 34 end 35end 36end%function 37 38function [Si,SiT] = bootstrap_h(fM1,fM2,fN,stg,scores_names) 39for n = 1:size(scores_names,2) 40 eval("fM1h=fM1." + scores_names(n) + ";"); 41 eval("fM2h=fM2." + scores_names(n) + ";"); 42 eval("fNh=fN." + scores_names(n) + ";"); 43 44 [Sih,SiTh] = calcSobolSaltelli(fM1h,fM2h,fNh,stg); 45 46 eval("Si." + scores_names(n) + "=Sih;"); 47 eval("SiT." + scores_names(n) + "=SiTh;"); 48end 49end 50 51function [Si,SiT] = bootstrap_hq(fM1,fM2,fN,stg,scores_names,j) 52for n = 1:size(scores_names,2) 53 eval("fM1h=fM1." + scores_names(n) + ";"); 54 eval("fM2h=fM2." + scores_names(n) + ";"); 55 eval("fNh=fN." + scores_names(n) + ";"); 56 57 rng(j*stg.rseed) 58 I=ceil(rand(size(fM1h,1),1)*size(fM1h,1)); 59 fM1q = fM1h(I,:); 60 fM2q = fM2h(I,:); 61 fNq = fNh(I,:,:); 62 63 [Sih,SiTh] = calcSobolSaltelli(fM1q,fM2q,fNq,stg); 64 eval("Si." + scores_names(n) + "=Sih;"); 65 eval("SiT." + scores_names(n) + "=SiTh;"); 66end 67end 68 69function [Si,SiT] = calcSobolSaltelli(fM1,fM2,fN,stg) 70%Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J. comp. neuroscience 27.3 (2009): 471.) 71 72[Nsamples,Nvars,Npars]=size(fN); 73 74if(stg.sasubmean) % Makes the model more stable 75 fM1 = fM1 - mean(fM1,1); 76 fM2 = fM2 - mean(fM2,1); 77 for i=1:Npars 78 fN(:,:,i) = fN(:,:,i) - mean(fN(:,:,i),1); 79 end 80end 81 82EY2 = mean(fM1.*fM2); % Valid definition (see Halnes et. al. Appendix) 83VY = sum(fM1.^2)/(Nsamples-1) - EY2; 84VYT = sum(fM2.^2)/(Nsamples-1) - EY2; 85 86Si = zeros(Nvars,Npars); 87SiT= zeros(Nvars,Npars); 88 89for i=1:Npars 90 Si(:,i) = (sum(fM1.*fN(:,:,i))/(Nsamples-1) - EY2)./VY; 91 SiT(:,i) = 1 - (sum(fM2.*fN(:,:,i))/(Nsamples-1) - EY2)./VYT; 92end 93end 94 95function rst = remove_sim_error(rst,stg) 96error=[]; 97error_helper=[]; 98 99for n = 1:size(rst.fM1.sd,1) 100 if max(rst.fM1.sd(n,:)) == stg.errorscore 101 error = [error,n]; 102 end 103 if max(rst.fM2.sd(n,:)) == stg.errorscore 104 error = [error,n]; 105 end 106 for m = 1:stg.parnum 107 if max(rst.fN.sd(n,:,m)) == stg.errorscore 108 error_helper = 1; 109 end 110 end 111 if error_helper == 1 112 error = [error,n]; 113 end 114 error_helper = 0; 115end 116 117error = unique(error); 118 119if ~isempty(error) 120 for n = size(error,2):-1:1 121 rst.fM1.sd(error(n),:) = []; 122 rst.fM2.sd(error(n),:) = []; 123 rst.fN.sd(error(n),:,:) = []; 124 125 rst.fM1.se(error(n),:) = []; 126 rst.fM2.se(error(n),:) = []; 127 rst.fN.se(error(n),:,:) = []; 128 129 rst.fM1.st(error(n),:) = []; 130 rst.fM2.st(error(n),:) = []; 131 rst.fN.st(error(n),:,:) = []; 132 133 rst.fM1.xfinal(error(n),:) = []; 134 rst.fM2.xfinal(error(n),:) = []; 135 rst.fN.xfinal(error(n),:,:) = []; 136 end 137end 138end
Takes the matrices fM1, fM2, and fN and calculates sensitivity indexes. It calculates indexes based on the following oputputs of the f_sim_score function:
Code modified from the Geir Halnes et al. 2009 paper.
References
Plots
f_plot
Code
1function f_plot(rst,stg,mmf) 2 3% Inform the user that the plots are being done 4disp("Plotting ...") 5 6data_model = mmf.model.data.data_model; 7 8% Import the data on the first run 9load(data_model,'Data','sbtab') 10 11% Generate figure with Scores 12if isfield(rst,'diag') 13 f_plot_scores(rst.diag,stg,sbtab) 14end 15 16% Generate figure with Inputs 17if isfield(rst,'diag') 18 f_plot_inputs(rst.diag,stg,sbtab) 19end 20 21% Generate figure with Outputs 22if isfield(rst,'diag') 23 f_plot_outputs(rst.diag,stg,sbtab,Data,mmf) 24end 25 26% Generate figure with input and Output of all experiments 27if isfield(rst,'diag') 28 f_plot_in_out(rst.diag,stg,sbtab,Data) 29end 30 31% Generate figure with optimization results 32if isfield(rst,'opt') 33 f_plot_opt(rst,stg) 34end 35 36% Generate figures for Sensitivity Analysis 37if isfield(rst,'gsa') 38 f_plot_gsa_sensitivities(rst.gsa,stg,sbtab); 39end 40 41% Generate figure for Profile Likelihood Analysis 42if isfield(rst,'PLA') 43 f_plot_PL(rst,stg,mmf) 44end 45 46end
The function that calls all the custom plot functions when appropriate Plots diagnosis that are important to understand if everything is working as it was supposed, it , expected outputs, observed outputs and scores for the models and conditions specified.
Outputs
Figure Scores
Total scores and scores per dataset given the parameters specified in stg.pa
Code Figure Scores
1function f_plot_scores(rst,stg,sbtab) 2set(0,'defaultTextFontName', 'Helvetica') 3set(0,'defaultAxesFontName', 'Helvetica') 4 5% Generates a figure with Scores 6 7% Inform the user that fig1 is being ploted 8disp("Plotting Scores") 9 10%Closes previous instances of the figure and generates a new one 11figHandles = findobj('type', 'figure', 'name', 'Scores'); 12close(figHandles); 13figure('WindowStyle', 'docked','Name','Scores','NumberTitle', 'off'); 14 15% Generate top plot of figure 1 16subplot(4,1,1) 17 18% Plot the total scores of each parameter array to test 19scatter(stg.pat,[rst(stg.pat).st],20,'filled') 20ylabel('Total Score ($s_t$)') 21set(gca,'xtick',[]) 22set(gca,'FontSize',10,'Fontweight','bold') 23 24% Choose correct title according to settings 25if stg.useLog == 1 26 title("Sum of the Log base 10 of the Score of each Experimental Output") 27elseif stg.useLog == 2 28 title("Sum of the Log base 10 of the Score of each Experiment") 29elseif stg.useLog == 3 30 title("Log base 10 of sum of the Score of all Experiments") 31elseif stg.useLog == 4 || stg.useLog == 0 32 title("Sum of the Score of all Experiments") 33end 34 35% Choose the bounds for the x axis so it aligns with the bottom plot 36xlim([min(stg.pat)-0.5 max(stg.pat)+0.5]) 37 38% Generate bottom plot of figure 1 39subplot(4,1,[2,3,4]) 40 41% Generate labels for left of heatmap (Experiment number Dataset number) 42label = []; 43 44% Iterate over the number of experiments 45for n = stg.exprun 46 47 % Iterate over the number of datasets in each experiment 48 for j = 1:size(sbtab.datasets(n).output,2) 49 50 label{size(label,2)+1} = {strrep("E" + (n-1) + " " + ... 51 string(sbtab.datasets(n).output_name{j}),"_","\_")}; 52 end 53end 54% Choose wether to use the score of each dataset or its log base 10 55% according to settings 56% if stg.useLog == 1 || stg.useLog == 4 57heatline = []; 58 59% Iterate over the number of parameter arrays to test 60for k = stg.pat 61 heatpoint{k} = []; 62 63 % Iterate over the number of experiments 64 for n = stg.exprun 65 66 % Get the score of each dataset 67 heatpoint{k} = [heatpoint{k};rst(k).sd(:,n)]; 68 69 end 70 % Combine heatpoints in order to correctly display heatmap 71 heatline = [heatline,heatpoint{k}]; 72end 73heatline( ~any(heatline,2), : ) = []; 74 75% Plot the heatmap 76h = heatmap(heatline,'Colormap',turbo,'YDisplayLabels',label,... 77 'GridVisible','off','FontSize',10); 78h.CellLabelFormat = '%.2e'; 79 80h.Title = "Score of each Experimental Output ($s_e$)"; 81h.XLabel = 'Parameter arrays ($\theta$)'; 82h.YLabel = 'Experimental Outputs'; 83 84h.NodeChildren(3).XAxis.Label.Interpreter = 'latex'; 85h.NodeChildren(3).YAxis.Label.Interpreter = 'latex'; 86% h.NodeChildren(3).ZAxis.Label.Interpreter = 'latex'; 87h.NodeChildren(3).Title.Interpreter = 'latex'; 88h.NodeChildren(3).TickLabelInterpreter = 'latex'; 89h.NodeChildren(2).TickLabelInterpreter = 'latex'; 90% h.NodeChildren(1).TickLabelInterpreter = 'latex'; 91% h.NodeChildren(1).Label.Interpreter = 'Latex'; 92% h.NodeChildren(2).Label.Interpreter = 'Latex'; 93 94end
Figure Inputs
Checks inputs to the model
Code Figure Inputs
1function f_plot_inputs(rst,stg,sbtab) 2% Generates a figure with Inputs, one subplot per experiment 3 4% Inform the user that fig2 is being ploted 5disp("Plotting Inputs") 6 7plot_n = 1; 8fig_n = 0; 9layout = []; 10% Iterate over the number of experiments 11for n = stg.exprun 12 13 % Generate the right amount of figures for all plots and calculates 14 % proper subploting position 15 16 [fig_n,layout] = f_get_subplot(size(stg.exprun,2),plot_n,fig_n,"Inputs",layout); 17 nexttile(layout); 18 19% fig_n = f_get_subplot(size(stg.exprun,2),plot_n,fig_n,"Inputs"); 20 21if mod(plot_n,24) == 1 22% Lgnd = legend('show','Orientation','Horizontal'); 23% Lgnd.Position(1) = 0; 24% Lgnd.Position(2) = 0.5; 25% Lgnd.Layout.Tile = 'North'; 26 xlabel(layout,"seconds", 'FontSize', 12,'Fontweight','bold','Interpreter','latex') 27% ylabel(layout,string(rst(m).simd{1,n}.DataInfo{end-... 28% size(sbtab.datasets(n).output,2)+j,1}.Units), 'FontSize', 12,'Fontweight','bold') 29% legend boxoff 30 31% ylabel(string(rst(m).simd{1,n}.DataInfo{end-... 32% size(sbtab.datasets(n).output,2)+j,1}.Units)) 33 end 34 35 36 plot_n = plot_n +1; 37 38 hold on 39 40 % Iterate over the number of inputs in each experiment 41 for j = 1:size(sbtab.datasets(n).input,2) 42 43 % Iterate over the number of parameter arrays to test 44 for m = stg.pat 45 46 % (Until a non broken simulation is found) 47 if rst(m).simd{1,n} ~= 0 48 49 % Plot the inputs to each experiment 50 plot(rst(m).simd{1,n}.Time,rst(m).simd{1,n}.Data(1:end,... 51 str2double(strrep(sbtab.datasets(n).input(j),'S',''))+1),'LineWidth',1.5) 52 53 % Get the correct label for each input of the experiment 54 labelfig2(j) = rst(m).simd{1,n}.DataNames(str2double(... 55 strrep(sbtab.datasets(n).input(j),'S',''))+1); 56 57 ylabel(layout,string(rst(m).simd{1,n}.DataInfo{... 58 str2double(strrep(sbtab.datasets(n).input(j),'S',''))+1,1}.Units),... 59 'FontSize', 12,'Fontweight','bold','Interpreter','latex') 60 61 break 62 end 63 end 64 end 65 66% xlabel('seconds') 67 % Add a legend to each plot 68 legend(labelfig2) 69 legend boxoff 70 clear labelfig2 71 72 ylim([0 inf]) 73 74 % Add a title to each plot 75 title("E"+(n-1)) 76 77 hold off 78end 79 80end
Figure Outputs
Expected outputs, observed outputs
Code Figure Outputs
1function f_plot_outputs(rst,stg,sbtab,Data,mmf) 2% Generates a figure with Outputs, one subplot per experimental output 3 4% Inform the user that fig3 is being ploted 5disp("Plotting Outputs") 6 7% Get the total number of outputs to set the total number of plots 8[plot_tn,~] = f_get_outputs(stg,sbtab); 9plot_n = 1; 10fig_n = 0; 11layout = []; 12% Iterate over the number of experiments 13for n = stg.exprun 14 15 % Iterate over the number of datasets in each experiment 16 for j = 1:size(sbtab.datasets(n).output,2) 17 18 % Generate the right amount of figures for all plots and calculates 19 % proper subploting position 20 21 [fig_n,layout] = f_get_subplot(plot_tn,plot_n,fig_n,"Outputs",layout); 22 nexttile(layout); 23 24% fig_n = f_get_subplot(plot_tn,plot_n,fig_n,"Outputs"); 25 26 % Add a legend to the figure 27 if mod(plot_n,24) == 1 28 Lgnd = legend('show','Orientation','Horizontal'); 29% Lgnd.Position(1) = 0; 30% Lgnd.Position(2) = 0.5; 31 Lgnd.Layout.Tile = 'North'; 32 xlabel(layout,"seconds", 'FontSize', 12,'Fontweight','bold','Interpreter','latex') 33% ylabel(layout,string(rst(m).simd{1,n}.DataInfo{end-... 34% size(sbtab.datasets(n).output,2)+j,1}.Units), 'FontSize', 12,'Fontweight','bold') 35 legend boxoff 36 37% ylabel(string(rst(m).simd{1,n}.DataInfo{end-... 38% size(sbtab.datasets(n).output,2)+j,1}.Units)) 39 end 40 41 plot_n = plot_n + 1; 42 43 hold on 44 45 % Iterate over the number of parameter arrays to test 46 for m = stg.pat 47 % (Until a non broken simulation is found) 48 if rst(m).simd{1,n} ~= 0 49 50 time = rst(m).simd{1,n}.Time; 51 data = Data(n).Experiment.x(:,j); 52 53 data_SD = Data(n).Experiment.x_SD(:,j); 54 55 % Plot the outputs to each dataset (new subplots) as they 56 % are given in the data provided in sbtab 57% scatter(time,data,'filled','k',... 58% 'DisplayName','data') 59 60 errorbar(time,data,data_SD,'ok','LineWidth',0.5,'MarkerSize',1,'DisplayName',"test"); 61 62 break 63 end 64 end 65 66 % Iterate over the number of parameter arrays to test 67 for m = stg.pat 68 69 % Plot only if the simulation was successful 70 if rst(m).simd{1,n} ~= 0 71 72 time = rst(m).simd{1,n}.Time; 73 [sim_results,~] = f_normalize(rst(m),stg,n,j,mmf); 74 if stg.simdetail 75 time_detailed = rst(m).simd{1,n+2*stg.expn}.Time; 76 [~,sim_results_detailed]= f_normalize(rst(m),stg,n,j,mmf); 77 end 78 79 % Plot the outputs to each dataset (new subplots) and 80 % parameter array to test that are simulated using 81 % Simbiology 82 if stg.simdetail 83 plot(time_detailed,... 84 sim_results_detailed,'DisplayName',... 85 string("$\theta_"+m + "$"),'LineWidth',1.5) 86 else 87 88 plot(time,... 89 sim_results,'DisplayName',... 90 string("$\theta_"+m + "$"),'LineWidth',1.5) 91 end 92 93% ylabel(string(rst(m).simd{1,n}.DataInfo{end-... 94% size(sbtab.datasets(n).output,2)+j,1}.Units)) 95 96 ylabel(layout,string(rst(m).simd{1,n}.DataInfo{end-... 97 size(sbtab.datasets(n).output,2)+j,1}.Units), 'FontSize', 12,'Fontweight','bold','Interpreter','latex') 98 end 99 end 100 101 hold off 102 103% xlabel('seconds') 104 105 if stg.simdetail 106 ylim([min([0,min(sim_results_detailed),min(sim_results),min(data-data_SD),min(data)]) inf]) 107 else 108 ylim([min([0,min(sim_results),min(data-data_SD),min(data)]) inf]) 109 end 110 111 % Choose correct title according to settings 112 if stg.plotoln == 1 113 title("E" + (n-1) + " " +... 114 strrep(string(sbtab.datasets(n).output_name{1,j}),'_','\_')) 115 else 116 title("E" + (n-1) + " " +... 117 string(sbtab.datasets(n).output{1,j})) 118 end 119 120 % Choose number of decimal places for y axis 121 ytickformat('%.2g') 122 end 123end 124end
Figure Input and Outputs per experiment
Combined figure of the inputs and outputs for each experiment, on the left side we have the inputs of the experiment and on the right side the outputs
Code Figure Input and Outputs
1function f_plot_in_out(rst,stg,sbtab,Data) 2% Generates a figure with input and Output of all experiments on the left 3% side it plots the inputs of the experiment and on the right side it plots 4% the outputs 5 6for n = stg.exprun 7 8 helper = 1; 9 f_plot_in_out_left(rst,stg,sbtab,helper,... 10 size(sbtab.datasets(n).output,2) > 4) 11 12 for j = 1:size(sbtab.datasets(n).output,2) 13 14 if j/4 > helper 15 helper = helper +1; 16 f_plot_in_out_left(rst,stg,sbtab,helper,... 17 size(sbtab.datasets(n).output,2) > 4) 18 end 19 20 if size(sbtab.datasets(n).output,2) == 1 21 subplot(1,2,j+ceil(j/(2/2))*1) 22 elseif size(sbtab.datasets(n).output,2) == 2 23 subplot(2,2,j+ceil(j/(2/2))*1) 24 elseif size(sbtab.datasets(n).output,2) > 2 &&... 25 size(sbtab.datasets(n).output,2) <= 4 26 subplot(2,4,j+ceil(j/(4/2))*2) 27 elseif size(sbtab.datasets(n).output,2) > 4 28 subplot(2,4,j+ceil(j/(4/2))*2-helper*8+8) 29 end 30 31 hold on 32 33 % Iterate over the number of parameter arrays to test 34 for m = stg.pat 35 % (Until a non broken simulation is found) 36 if rst(m).simd{1,n} ~= 0 37 38 % Plot the outputs to each dataset (new subplots) as they 39 % are given in the data provided in sbtab 40 41 time = rst(m).simd{1,n}.Time; 42 data = Data(n).Experiment.x(:,j); 43 data_SD = Data(n).Experiment.x_SD(:,j); 44 45 scatter(time,data,'filled','k',... 46 'DisplayName','data') 47 48 errorbar(time,data,data_SD, 'vertical', 'k', 'LineStyle', 'none','LineWidth',1); 49 break 50 end 51 end 52 % Iterate over the number of parameter arrays to test 53 for m = stg.pat 54 55 % Plot only if the simulation was successful 56 if rst(m).simd{1,n} ~= 0 57 58 time = rst(m).simd{1,n}.Time; 59 [sim_results] = f_normalize(rst(m),stg,n,j); 60 61 62 if stg.simdetail 63 time_detailed = rst(m).simd{1,n+2*stg.expn}.Time; 64 [~,sim_results_detailed]= f_normalize(rst(m),stg,n,j); 65 end 66 67 % Plot the outputs to each dataset (new subplots) and 68 % parameter array to test that are simulated using 69 % Simbiology 70 if stg.simdetail 71 plot(time_detailed,... 72 sim_results_detailed,'DisplayName',... 73 string("Parameter set "+m),'LineWidth',1.5) 74 else 75 plot(time,... 76 sim_results,... 77 'DisplayName',string("Parameter set "+m),... 78 'LineWidth',1.5) 79 end 80 81 ylabel(string(rst(m).simd{1,n}.DataInfo{end-... 82 size(sbtab.datasets(n).output,2)+j,1}.Units),... 83 'FontSize', 12,'Fontweight','bold') 84 end 85 end 86 87 hold off 88 89 set(gca,'FontSize',12,'Fontweight','bold') 90 91 if stg.simdetail 92 ylim([min([0,min(sim_results_detailed),min(sim_results),min(data-data_SD),min(data)]) inf]) 93 else 94 ylim([min([0,min(sim_results),min(data-data_SD),min(data)]) inf]) 95 end 96 97 xlabel('seconds','FontSize', 12,'Fontweight','bold') 98 99 % Choose correct title according to settings 100 if stg.plotoln == 1 101 title(strrep(string(sbtab.datasets(n).output_name{1,j}),'_',... 102 '\_'),'FontSize', 16,'Fontweight','bold') 103 else 104 title(string(sbtab.datasets(n).output{1,j}),'FontSize', 16,... 105 'Fontweight','bold') 106 end 107 108 % Choose number of decimal places for y axis 109 ytickformat('%.2g') 110 end 111end 112 113 function f_plot_in_out_left(rst,stg,sbtab,helper,reuse) 114 if reuse 115 figHandles = findobj('type', 'figure', 'name', "E " + (n-1) +... 116 " " + helper); 117 close(figHandles); 118 figure('WindowStyle', 'docked','Name', "E " + (n-1)+ " " +... 119 helper,'NumberTitle', 'off'); 120 sgtitle( "Experiment " + (n-1) + " " + helper + " (E " +... 121 (n-1) + " " + helper +")",'FontSize', 28); 122 else 123 figHandles = findobj('type', 'figure', 'name', "E " + (n-1)); 124 close(figHandles); 125 figure('WindowStyle', 'docked','Name', "E " + (n-1),... 126 'NumberTitle', 'off'); 127 sgtitle( "Experiment " + (n-1) + " (E " + (n-1) +... 128 ")",'FontSize', 28); 129 end 130 131 subplot(2,4,[1,2,5,6]) 132 133 hold on 134 for o = 1:size(sbtab.datasets(n).input,2) 135 for p = stg.pat 136 137 % (Until a non broken simulation is found) 138 if rst(p).simd{1,n} ~= 0 139 % Plot the inputs to each experiment 140 plot(rst(p).simd{1,n}.Time,rst(p).simd{1,n}.Data(1:end,... 141 str2double(strrep(sbtab.datasets(n).input(o),'S','')... 142 )+1),'LineWidth',1.5) 143 144 % Get the correct label for each input of the 145 % experiment 146 labelfig2(o) = rst(p).simd{1,n}.DataNames(str2double(... 147 strrep(sbtab.datasets(n).input(o),'S',''))+1); 148 149 ylabel(string(rst(p).simd{1,n}.DataInfo{... 150 str2double(strrep(sbtab.datasets(n).input(o),'S',''))+1,1}.Units),... 151 'FontSize', 12,'Fontweight','bold') 152 153 154 break 155 end 156 end 157 158 end 159 160 set(gca,'FontSize',12,'Fontweight','bold') 161 162 xlabel('seconds','FontSize', 12,'Fontweight','bold') 163 % Add a legend to each plot 164 legend(labelfig2,'FontSize', 16,'Fontweight','bold') 165 legend boxoff 166 clear labelfig2 167 168 ylim([0 inf]) 169 170 % Add a title to each plot 171 title("Inputs",'FontSize', 18,'Fontweight','bold') 172 173 hold off 174 end 175end
Figure Sensitivity Analysis \(S_{i}\)
Figure Sensitivity Analysis \(S_{Ti}\)
Code figures SA
1function f_plot_gsa_sensitivities(rst,stg,sbtab) 2% Generates figures for Sensitivity Analysis 3 4% Get the total number of outputs 5[~,outputNames.sd] = f_get_outputs(stg,sbtab); 6 7for n = 1:size(outputNames.sd,2) 8 outputNames.sd{n}{:} = strrep(outputNames.sd{n}{:},"_","\_"); 9end 10for n = stg.exprun 11 outputNames.se{n} = "E " + string(n-1); 12end 13 14outputNames.xfinal = outputNames.sd; 15 16parNames = cell(1,stg.parnum); 17parNames2 = cell(1,stg.parnum); 18 19for n = 1:stg.parnum 20 parNames{n} = char("P" + find(stg.partest==n)); 21end 22 23for n = 1:size(parNames,2) 24 parNames2{n} = string(parNames{n}(1,:)); 25 for m = 2:size(parNames{n},1) 26 parNames2{n} = string(parNames2{n}) + ", " +... 27 string(parNames{n}(m,:)); 28 end 29end 30 31% Bootstrapping quartile mean of first order Sensitivity index for score 32% per Experimental Output 33f_generate_plot(rst,stg,outputNames,parNames2,... 34 "Si seo bm",... 35 "First order Sensitivities calculated using the Score of each Experimental Output (Bootstrapping Mean)",... 36 "outputNames.sd",... 37 "transpose(reshape(mean(rst.SiQ.sd(:,:,1:stg.parnum)),[size(rst.SiQ.sd,2),size(rst.SiQ.sd,3)]))") 38 39% Bootstrapping quartile mean of total order Sensitivity index for score 40% per Experimental Output 41f_generate_plot(rst,stg,outputNames,parNames2,"SiT seo bm",... 42 "Total order Sensitivities calculated using the Score of each Experimental Output (Bootstrapping Mean)",... 43 "outputNames.sd",... 44 "transpose(reshape(mean(rst.SiTQ.sd(:,:,1:stg.parnum)),[size(rst.SiQ.sd,2),size(rst.SiQ.sd,3)]))") 45 46% Bootstrapping quartile mean of first order Sensitivity index for score 47% per Experiments 48f_generate_plot(rst,stg,outputNames,parNames2,"Si se bm",... 49 "First order Sensitivities calculated using the Score of each Experiment(Bootstrapping Mean)",... 50 "outputNames.se",... 51 "transpose(reshape(mean(rst.SiQ.se(:,:,1:stg.parnum)),[size(rst.SiQ.se,2),size(rst.SiQ.se,3)]))") 52 53% Bootstrapping quartile mean of total order Sensitivity index for score 54% per Experiments 55f_generate_plot(rst,stg,outputNames,parNames2,"SiT se bm",... 56 "Total order Sensitivities calculated using the Score of each Experiment (Bootstrapping Mean)",... 57 "outputNames.se",... 58 "transpose(reshape(mean(rst.SiTQ.se(:,:,1:stg.parnum)),[size(rst.SiQ.se,2),size(rst.SiQ.se,3)]))") 59 60% Bootstrapping quartile mean of first order Sensitivity index for the 61% final points of the simulations for the output beeing measured 62f_generate_plot(rst,stg,outputNames,parNames2,"Si xfinal bm",... 63 "First order Sensitivities calculated using the final value of each Experimental Output (Bootstrapping Mean)",... 64 "outputNames.xfinal",... 65 "transpose(reshape(mean(rst.SiQ.xfinal(:,:,1:stg.parnum)),[size(rst.SiQ.xfinal,2),size(rst.SiQ.xfinal,3)]))") 66 67% Bootstrapping quartile mean of total order Sensitivity index for the 68% final points of the simulations for the output beeing measured 69f_generate_plot(rst,stg,outputNames,parNames2,"SiT xfinal bm",... 70 "Total order Sensitivities calculated using the final value of each Experimental Output (Bootstrapping Mean)",... 71 "outputNames.xfinal",... 72 "transpose(reshape(mean(rst.SiTQ.xfinal(:,:,1:stg.parnum)),[size(rst.SiQ.xfinal,2),size(rst.SiQ.xfinal,3)]))") 73 74figHandles = findobj('type', 'figure', 'name', 'Si,SiT'); 75close(figHandles); 76figure('WindowStyle', 'docked','Name','Si,SiT', 'NumberTitle', 'off'); 77 78for n = 1:size(parNames2,2) 79 a{n} = char(parNames2{n}); 80end 81 82a = categorical(a,a); 83 84bar(a,[transpose(rst.Si.st(:,1:stg.parnum)),... 85 transpose(rst.SiT.st(:,1:stg.parnum))]) 86xlabel('Parameters'); 87ylabel('Sensitivity'); 88title('Sensitivities calculated using the sum of the Score of all Experiments'); 89legend({'SI','SIT'}); 90legend boxoff 91 92figHandles = findobj('type', 'figure', 'name', 'Si,SiT b'); 93close(figHandles); 94figure('WindowStyle', 'docked','Name','Si,SiT b', 'NumberTitle', 'off'); 95 96T = []; 97 98for n = 1:size(a,2) 99 for m = 1:size(rst.SiQ.st(:,n),1) 100 T = [T;table(rst.SiQ.st(m,n),a(n),"Si")]; 101 end 102end 103for n = 1:size(a,2) 104 for m = 1:size(rst.SiTQ.st(:,n),1) 105 T = [T;table(rst.SiTQ.st(m,n),a(n),"SiT")]; 106 end 107end 108 109boxchart(T.Var2,T.Var1,'GroupByColor',T.Var3,'MarkerStyle','.','JitterOutliers','on') 110xlabel('Parameters'); 111ylabel('Sensitivity'); 112title('Sensitivities calculated using the sum of the Score of all Experiments (Bootstrapping)'); 113legend({'Si','SiT'},'Location','best'); 114legend boxoff 115end 116 117function f_generate_plot(rst,stg,outputNames,parNames2,name,title,... 118 helprer1,helprer2) 119 120eval("figHandles = findobj('type', 'figure', 'name', '" + name + "');") 121close(figHandles); 122eval("figure('WindowStyle', 'docked','Name','" + name +... 123 "','NumberTitle', 'off');") 124 125heatmap_fixer = eval(helprer1); 126heatmap_fixer=heatmap_fixer(~cellfun('isempty',heatmap_fixer)); 127 128heatmap_fixer2 = eval(helprer2); 129heatmap_fixer2 = heatmap_fixer2(:,all(~isnan(heatmap_fixer2))); 130 131h = heatmap(heatmap_fixer,parNames2,heatmap_fixer2,'Colormap',turbo,... 132 'ColorLimits',[0 1],'GridVisible','off'); 133h.CellLabelFormat = '%.2f'; 134 135eval(" h.Title = """ + title + """;") 136h.XLabel = 'Outputs'; 137h.YLabel = 'Parameters'; 138end
Calls
Loads - data.mat
f_get_subplot
Code
1function [fig_n,layout] = f_get_subplot(plot_tn,plot_n,fig_n,fig_name,layout) 2 3% size_x = 4; 4% size_y = 6; 5size_t = 24; 6% ratio_1 = 3; 7% ratio_2 = 4; 8 9 10size_x = [1,1,1,2,3,3,4,4,3,4,4,4,5,5,5,4,5,5,5,5,6,6,6,6]; 11size_y = [1,2,3,2,2,2,2,2,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4]; 12 13% If the amount of plots is bigger thatn the maximum amount of plots per 14% figure subdivide the plots to more than one figure 15% if plot_tn > 24 16% % Generate a new figure for the first plot and each time the number of 17% % plots is greater than figure number divided by max plot number per 18% % figure 19% if mod(plot_n-1,24) == 0 20% 21% fig_n = fig_n + 1; 22% 23% %Close previous instances of the figure and generates a new one 24% figHandles = findobj('type', 'figure', 'name', fig_name + " " + fig_n); 25% close(figHandles); 26% figure('WindowStyle', 'docked','Name', fig_name + " " + fig_n,'NumberTitle', 'off'); 27% sgtitle(fig_name + " " + fig_n); 28% layout = tiledlayout(ceil(sqrt(24/6)*2),ceil(sqrt(24/6)*3),'Padding','none','TileSpacing','compact'); 29% end 30% 31% % Get the correct subploting position for each plot 32% % if plot_tn/24 < fig_n 33% % subplot(ceil(sqrt((plot_tn-(fig_n-1)*24)/6)*2),ceil(sqrt((plot_tn-(fig_n-1)*24)/6)*3),plot_n-(fig_n-1)*24) 34% % else 35% % subplot(ceil(sqrt(24/6)*2),ceil(sqrt(24/6)*3),plot_n-(fig_n-1)*24) 36% % end 37% else 38 39 % Generate a new figure for the first plot 40 if mod(plot_n-1,size_t) == 0 41 42 %Close previous instances of the figure and generates a new one 43 figHandles = findobj('type', 'figure', 'name', fig_name); 44 close(figHandles); 45 figure('WindowStyle', 'docked','Name', fig_name, 'NumberTitle', 'off'); 46 sgtitle(fig_name); 47% layout = tiledlayout(... 48% ceil(sqrt((plot_tn-fig_n*size_t)/(ratio_1*ratio_2))*ratio_1),... 49% ceil(sqrt((plot_tn-fig_n*size_t)/(ratio_1*ratio_2))*ratio_2),... 50% 'Padding','none','TileSpacing','compact'); 51 layout = tiledlayout(... 52 size_x(plot_tn-(floor(plot_tn/size_t)*size_t)),... 53 size_y(plot_tn-(floor(plot_tn/size_t)*size_t)),... 54 'Padding','none','TileSpacing','compact'); 55 56% title(layout,fig_name, 'FontSize', 16,'Fontweight','bold') 57 fig_n = fig_n + 1; 58 59 end 60 61 % Get the correct subploting position for each plot 62% subplot(ceil(sqrt(plot_tn/6)*2),ceil(sqrt(plot_tn/6)*3),plot_n) 63 64% end 65end
Inputs
Outputs
Calls
Loads
General purpose
f_save_analysis
Code
1function f_save_analysis(stg,sb,rst,mmf) 2 3Results_Folder = mmf.model.results.main; 4Analysis_folder = mmf.model.results.analysis.main; 5Analysis_date_folder = mmf.model.results.analysis.date.main; 6 7[~,~] = mkdir(Results_Folder); 8[~,~] = mkdir(Analysis_folder); 9[~,~] = mkdir(Analysis_date_folder); 10addpath(Analysis_date_folder) 11 12save (Analysis_date_folder + "Analysis.mat",'stg','sb','rst'); 13end
Inputs
Outputs
Calls
Loads
f_save_plots
Code
1function f_save_plots(mmf) 2 3Analysis_date_folder = mmf.model.results.analysis.date.main; 4 5FigList = findobj(allchild(0), 'flat', 'Type', 'figure'); 6 7[~,~] = mkdir(Analysis_date_folder); 8 9savefig(FigList(end:-1:1),... 10 Analysis_date_folder + "All_figures.fig"); 11 12for iFig = 1:length(FigList) 13 FigHandle = FigList(iFig); 14 FigName = get(FigHandle, 'Name'); 15 16 saveas(FigHandle, Analysis_date_folder + FigName + ".png") 17end 18end
Inputs
Outputs
Calls
Loads
f_get_outputs
Code
1function [nOutputs,outputNames] = f_get_outputs(stg,sbtab) 2 3persistent n_out 4persistent out_name 5 6if isempty(n_out) 7 n_out = 0; 8 out_name = []; 9 for n = stg.exprun 10 for j = 1:size(sbtab.datasets(n).output,2) 11 n_out = n_out + 1; 12 out_name{n_out} = {"E" + (n-1) + " " + string(sbtab.datasets(n).output{1,j})}; 13 end 14 end 15end 16 17nOutputs = n_out; 18outputNames = out_name; 19end
Inputs
Outputs
Calls
Loads
Settings file
Default settings code
1function [stg] = default_settings() 2 3%% Import 4 5% True or false to decide whether to run import functions (Import) 6stg.import = true; 7 8 9% Name of the excel file with the sbtab (SBtab excel name) 10stg.sbtab_excel_name = "SBTAB.xlsx"; 11 12% Name of the model (Name) 13stg.name = "model_name"; 14 15% Name of the default model compartment (Compartment name) 16stg.cname = "Compartment"; 17 18% Name of the sbtab saved in .mat format (SBtab name) 19stg.sbtab_name = "SBtab_" + stg.name; 20 21%% Analysis 22 23% Experiments to run (example experiment 1 to 3 and experimet 6) 24stg.exprun = [1:3,6]; 25 26% Choice between 0,1,2 and 3,4 to choose the scoring method (check 27% documentation) (Use logarithm) 28stg.useLog = 4; 29 30% True or false to decide whether to use multicore everywhere it is 31% available (Optimization Multicore) 32stg.optmc = true; 33 34% Choice of ramdom seed (Ramdom seed) 35stg.rseed = 1; 36 37% True or false to decide whether to use display simulation diagnostics in 38% the console (Simulation Console) 39stg.simcsl = false; 40 41% True or false to decide whether to display optimization results on 42% console (Optimization console) 43stg.optcsl = false; 44 45% True or false to decide whether to save results (Save results) 46stg.save_results = true; 47 48% True or false to decide whether to run detailed simulation for plots 49stg.simdetail = true; 50 51%% Simulation 52 53% Maximum time for each individual function to run in seconds (Maximum 54% time) 55stg.maxt = 10; 56 57% Equilibration time (Equilibration time) 58stg.eqt = 50000; 59 60% True or false to decide whether to do Dimensional Analysis (Dimensional 61% Analysis) 62stg.dimenanal = false; 63 64% True or false to decide whether to do Unit conversion (Unit conversion) 65stg.UnitConversion = false; 66 67% True or false to decide whether to do Absolute Tolerance Scaling 68% (Absolute Tolerance Scaling) 69stg.abstolscale = true; 70 71% Value of Relative tolerance (Relative tolerance) 72stg.reltol = 1.0E-4; 73 74% Value of Absolute tolerance (Absolute tolerance) 75stg.abstol = 1.0E-4; 76 77% Time units for simulation (Simulation time) 78stg.simtime = "second"; 79 80% True or false to decide whether to run sbioaccelerate (after changing 81% this value you need to run "clear functions" to see an effect) 82% (sbioaccelerate) 83stg.sbioacc = true; 84 85% (Absolute tolerance step size for equilibration) 86stg.abstolstepsize_eq = []; 87 88% Max step size in the simulation (if empty matlab decides whats best) 89% (Maximum step) 90stg.maxstep = [10]; 91 92% Max step size in the equilibration (if empty matlab decides whats best) 93% (Maximum step) 94stg.maxstepeq = [2]; 95 96% Max step size in the detailed plots (if empty matlab decides whats best) 97% (Maximum step) 98stg.maxstepdetail = [2]; 99 100% Default score when there is a simulation error, this is needed to keep 101% the optimizations working. 102% (error score) 103stg.errorscore = 10^5; 104%% Model 105 106% Number of parameters to optimize (Parameter number) 107stg.parnum = 5; 108 109original_parameter_set = zeros(1,10); 110 111% Array with the lower bound of all parameters (Lower bound) 112stg.lb = original_parameter_set-5; 113 114% Array with the upper bound of all parameters (Upper bound) 115stg.ub = original_parameter_set+5; 116 117%% Diagnostics 118 119% Choice of what parameters in the array to test, the indices correspond to 120% the parameters in the model and the numbers correspond to the parameters 121% in the optimization array, usually not all parameters are optimized so 122% there needs to be a match between one and the other. (Parameters to test) 123% In this example there are ten parameters in this imaginary model and we 124% are only interested in parameter 2,4,8,9, and 10. Note that stg.parnum is 125% five because of this and not ten 126stg.partest(:,1) = [0,1,0,2,0,0,0,3,4,5]; 127 128% (Parameter array to test) 129stg.pat = [1:2]; 130 131% All the parameter arrays, in this case there is only one (parameters here 132% are in log10 space)(Parameter arrays) 133stg.pa(1,:) = [1,1,1,1,1]; 134stg.pa(1,:) = [1,0,1,2,1]; 135 136% Best parameter array found so far for the model (Best parameter array) 137stg.bestpa = stg.pa(1,:); 138 139%% Plots 140 141% True or false to decide whether to plot results (Plots) 142stg.plot = true; 143 144% True or false to decide whether to use long names in the title of the 145% outputs plots in f_plot_outputs.m (Plot outputs long names) 146stg.plotoln = true; 147 148%% Sensitivity analysis 149 150% Number of samples to use in SA (Sensitivity analysis number of samples) 151stg.sansamples = 100; 152 153% True or false to decide whether to subtract the mean before calculating 154% SI and SIT (Sensitivity analysis subtract mean) 155stg.sasubmean = true; 156 157% Choose the way you want to obtain the samples of the parameters for 158% performing the SA; 0 Log uniform distribution truncated at the parameter 159% bounds 1 Log normal distribution with mu as the best value for a 160% parameter and sigma as stg.sasamplesigma truncated at the parameter 161% bounds 2 same as 1 without truncation 3 Log normal distribution centered 162% at the mean of the parameter bounds and sigma as stg.sasamplesigma 163% truncated at the parameter bounds 4 same as 3 without truncation. 164% (Sensitivity analysis sampling mode) 165stg.sasamplemode = 2; 166 167% Sigma for creating the normal distribution of parameters to perform 168% sensitivity analysis (Sensitivity analysis sampling sigma) 169stg.sasamplesigma = 0.1; 170 171%% Optimization 172 173% Time for the optimization in seconds (fmincon does not respect this 174% time!!) (Optimization time) 175stg.optt = 60*5; 176 177% Population size for the algorithms that use populations (Population size) 178stg.popsize = 144; 179 180% optimization start method, choose between: 1 Random starting point or 181% group of starting points inside the bounds 2 Random starting point or 182% group of starting points near the best point (Optimization start method) 183stg.osm = 1; 184 185% Distance from best parameter array to be used in stg.osm method 2 186% (Distance from best parameter array) 187stg.dbs = 0.1; 188 189% True or false to decide whether to use Multistart (Multistart) 190stg.mst = false; 191 192% Multistart size 193stg.msts = 1; 194 195% True or false to decide whether to display Plots (Plots doesn't work if 196% using multicore) (Optimization plots) 197stg.optplots = true; 198 199% True or false to decide whether to run fmincon (no gradient so this 200% doesn't work very well, no max time!!) 201stg.fmincon = false; 202 203% Options for fmincon (fmincon options) 204stg.fm_options = optimoptions('fmincon',... 205 'Algorithm','interior-point',... 206 'MaxIterations',2,'OptimalityTolerance',0,... 207 'StepTolerance',1e-6,'FiniteDifferenceType','central',... 208 'MaxFunctionEvaluations',10000); 209 210% True or false to decide whether to run simulated annealing (Simulated 211% annealing) 212stg.sa = false; 213 214% Options for simulated annealing (Simulated annealing options) 215stg.sa_options = optimoptions(@simulannealbnd, ... 216 'MaxTime',stg.optt,... 217 'ReannealInterval',40); 218 219% True or false to decide whether to run Pattern search (Pattern search) 220stg.psearch = false; 221 222% Options for Pattern search (Pattern search options) 223%stg.psearch_options = optimoptions(@patternsearch,... 224 %'MaxTime',stg.optt,'UseParallel',stg.optmc,... 225 %'UseCompletePoll',true,'UseCompleteSearch',true,... 226 %'MaxMeshSize',2,'MaxFunctionEvaluations',2000); 227 228% True or false to decide whether to run Genetic algorithm (Genetic 229% algorithm) 230stg.ga = false; 231 232% Options for Genetic algorithm (Genetic algorithm options) 233stg.ga_options = optimoptions(@ga,'MaxGenerations',200,... 234 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 235 'PopulationSize',stg.popsize,... 236 'MutationFcn','mutationadaptfeasible'); 237 238% True or false to decide whether to run Particle swarm (Particle swarm) 239stg.pswarm = false; 240 241% Options for Particle swarm (Particle swarm options) 242stg.pswarm_options = optimoptions('particleswarm',... 243 'MaxTime',stg.optt,'UseParallel',stg.optmc,'MaxIterations',200,... 244 'SwarmSize',stg.popsize); 245 246% True or false to decide whether to run Surrogate optimization (Surrogate 247% optimization) 248stg.sopt = false; 249 250% Options for Surrogate optimization (Surrogate optimization options) 251stg.sopt_options = optimoptions('surrogateopt',... 252 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 253 'MaxFunctionEvaluations',5000,... 254 'MinSampleDistance',0.2,'MinSurrogatePoints',32*2+1); 255endExample settings code
1function [stg] = Example_model() 2 3%% Import 4 5% True or false to decide whether to run import functions (Import) 6stg.import = true; 7 8% Name of the excel file with the sbtab (SBtab excel name) 9stg.sbtab_excel_name = "SBTAB example.xlsx"; 10 11% Name of the model (Name) 12stg.name = "Example"; 13 14% Name of the default model compartment (Compartment name) 15stg.cname = "Compartment"; 16 17%% Analysis 18 19% Experiments to run 20stg.exprun = [1,2]; 21% stg.exprun = [1,2,3]; 22 23% Choice between 0,1,2 and 3 to change either and how to apply log10 to the 24% scores (check documentation) (Use logarithm) 25stg.useLog = 0; 26 27% True or false to decide whether to use multicore everywhere it is 28% available (Optimization Multicore) 29stg.optmc = false; 30 31% Choice of ramdom seed (Ramdom seed) 32stg.rseed = 1; 33 34% True or false to decide whether to use display simulation diagnostics in 35% the console (Simulation Console) 36stg.simcsl = false; 37 38% True or false to decide whether to display optimization results on 39% console (Optimization console) 40stg.optcsl = true; 41 42% True or false to decide whether to display PLA results on console (PLA 43% console) 44stg.placsl = true; 45 46% True or false to decide whether to save results (Save results) 47stg.save_results = true; 48 49% True or false to decide whether to run detailed simulation for plots 50stg.simdetail = false; 51 52%% Simulation 53 54% Maximum time for each individual function to run in seconds (Maximum 55% time) 56stg.maxt = 2; 57 58% Equilibration time (Equilibration time) 59stg.eqt = 50000; 60 61% True or false to decide whether to do Dimensional Analysis (Dimensional 62% Analysis) 63stg.dimenanal = true; 64 65% True or false to decide whether to do Unit conversion (Unit conversion) 66stg.UnitConversion = true; 67 68% True or false to decide whether to do Absolute Tolerance Scaling 69% (Absolute Tolerance Scaling) 70stg.abstolscale = true; 71 72% Value of Relative tolerance (Relative tolerance) 73stg.reltol = 1.0E-4; 74 75% Value of Absolute tolerance (Absolute tolerance) 76stg.abstol = 1.0E-7; 77 78% Time units for simulation (Simulation time) 79stg.simtime = "second"; 80 81% True or false to decide whether to run sbioaccelerate (after changing 82% this value you need to run "clear functions" to see an effect) 83% (sbioaccelerate) 84stg.sbioacc = false; 85 86% Max step size in the simulation (if empty matlab decides whats best) 87% (Maximum step) 88stg.maxstep = []; 89 90% Max step size in the equilibration (if empty matlab decides whats best) 91% (Maximum step) 92stg.maxstepeq = []; 93 94% Max step size in the detailed plots (if empty matlab decides whats best) 95% (Maximum step) 96stg.maxstepdetail = [0.001]; 97 98% Default score when there is a simulation error, this is needed to keep 99% the optimizations working. (error score) 100stg.errorscore = 10^5; 101 102%% Model 103 104% Number of parameters to optimize (Parameter number) 105stg.parnum = 12; 106 107% Index for the parameters that have thermodynamic constrains (Termodiamic 108% Constrains Index) 109stg.tci = [8]; 110 111% Parameters to multiply to the first parameter (in Stg.partest to get to 112% the correct thermodynamic constrain formula) (Termodiamic Constrains 113% multipliers) 114stg.tcm([8],1) = [4]; 115stg.tcm([8],2) = [5]; 116stg.tcm([8],3) = [7]; 117 118% Parameters to divide to the first parameter (in Stg.partest to get to the 119% correct thermodynamic constrain formula) (Termodiamic Constrains 120% divisors) 121stg.tcd([8],1) = [1]; 122stg.tcd([8],2) = [3]; 123stg.tcd([8],3) = [6]; 124 125% Array with the lower bound of all parameters (Lower bound) 126stg.lb = zeros(1,stg.parnum)-4; 127 128% Array with the upper bound of all parameters (Upper bound) 129stg.ub = zeros(1,stg.parnum)+4; 130 131%% Diagnostics 132 133% Choice of what parameters in the array to test, the indices correspond to 134% the parameters in the model and the numbers correspond to the parameters 135% in the optimization array, usually not all parameters are optimized so 136% there needs to be a match between one and the other. (Parameters to test) 137stg.partest(:,1) = [1 ,2 ,3 ,4 ,5 ,6 ,7 ,2 ,8 ,9,... 138 10 ,11 ,12]; 139 140% (Parameter array to test) 141stg.pat = 1:3; 142 143% All the parameter arrays, in this case there is only one (Parameter 144% arrays) 145stg.pa(1,:) = [3.999,0.696,1.072,3.429,-0.751,-3.741,-0.569,0.831,3.068,0.921,-2.156,-1.970]; 146stg.pa(2,:) = stg.pa(1,:)-1; 147stg.pa(3,:) = stg.pa(1,:)+1; 148 149% Best parameter array found so far for the model (Best parameter array) 150stg.bestpa = stg.pa(1,:); 151 152%% Plots 153 154% True or false to decide whether to plot results (Plots) 155stg.plot = true; 156 157% True or false to decide whether to use long names in the title of the 158% outputs plots in f_plot_outputs.m (Plot outputs long names) 159stg.plotnames = true; 160 161%% Sensitivity analysis 162 163% Number of samples to use in SA (Sensitivity analysis number of samples) 164stg.sansamples = 1000; 165 166% True or false to decide whether to subtract the mean before calculating 167% SI and SIT (Sensitivity analysis subtract mean) 168stg.sasubmean = true; 169 170% Choose the way you want to obtain the samples of the parameters for 171% performing the SA; 0 Log uniform distribution truncated at the parameter 172% bounds 1 Log normal distribution with mu as the best value for a 173% parameter and sigma as stg.sasamplesigma truncated at the parameter 174% bounds 2 same as 1 without truncation 3 Log normal distribution centered 175% at the mean of the parameter bounds and sigma as stg.sasamplesigma 176% truncated at the parameter bounds 4 same as 3 without truncation. 177% (Sensitivity analysis sampling mode) 178stg.sasamplemode = 2; 179 180% Sigma for creating the normal distribution of parameters to perform 181% sensitivity analysis (Sensitivity analysis sampling sigma) 182stg.sasamplesigma = 0.1; 183 184stg.gsabootstrapsize = ceil(sqrt(stg.sansamples)); 185 186%% Profile Likelihood 187 188% Parameter(optimization array) that is being worked on in a specific 189% iteration of PL (if -1 no parameter is being worked in PL) 190% (Profile Likelihood Index) 191stg.PLind = -1; 192 193% Which parameters to do PL on, it should be all parameters but can also be 194% a subset for testing purposes 195% (Profile Likelihood parameters to Test) 196stg.pltest = (1:12); 197 198% How many points to do for each parameter in the PL 199% (Profile Likelihood Resolution) 200stg.plres = 80; 201 202% True or false to decide whether to do plots after calculating PL 203% (Profile Likelihood Plots) 204stg.plplot = true; 205 206% True or false to decide whether to run simulated annealing 207% (Profile Likelihood Simulated Annealing) 208stg.plsa = true; 209 210% Options for simulated annealing 211stg.plsao = optimoptions(@simulannealbnd,'Display','off', ... 212 'InitialTemperature',... 213 ones(1,stg.parnum)*1,'MaxTime',5,'ReannealInterval',40); 214 215% 0 or 1 to decide whether to run fmincon 216% (Profile Likelihood FMincon) 217stg.plfm = false; 218 219% Options for fmincon 220stg.plfmo = optimoptions('fmincon','Display','off',... 221 'Algorithm','interior-point',... 222 'MaxIterations',5,'OptimalityTolerance',0,... 223 'StepTolerance',1e-6,'FiniteDifferenceType','central'); 224 225%% Optimization 226 227% Time for the optimization in seconds (fmincon does not respect this 228% time!!) (Optimization time) 229stg.optt = 60*1; 230 231% Population size for the algorithms that use populations (Population size) 232stg.popsize = 10; 233 234% optimization start method, choose between: 1 Random starting point or 235% group of starting points inside the bounds 2 Random starting point or 236% group of starting points near the best point (Optimization start method) 237stg.osm = 1; 238 239% Distance from best parameter array to be used in stg.osm method 2 240% (Distance from best parameter array) 241stg.dbs = 0.1; 242 243% True or false to decide whether to use Multistart (Multistart) 244stg.mst = false; 245 246% Multistart size 247stg.msts = 1; 248 249% True or false to decide whether to display Plots (Plots doesn't work if 250% using multicore) (Optimization plots) 251stg.optplots = true; 252 253% True or false to decide whether to run fmincon (no gradient so this 254% doesn't work very well, no max time!!) 255stg.fmincon = false; 256 257% Options for fmincon (fmincon options) 258stg.fm_options = optimoptions('fmincon',... 259 'UseParallel',stg.optmc,... 260 'Algorithm','interior-point',... 261 'MaxIterations',2,'OptimalityTolerance',0,... 262 'StepTolerance',1e-6,'FiniteDifferenceType','central',... 263 'MaxFunctionEvaluations',10000); 264 265% True or false to decide whether to run simulated annealing (Simulated 266% annealing) 267stg.sa = false; 268 269% Options for simulated annealing (Simulated annealing options) 270stg.sa_options = optimoptions(@simulannealbnd, ... 271 'MaxTime',stg.optt,... 272 'InitialTemperature',... 273 ones(1,stg.parnum)*2,'ReannealInterval',40); 274 275% True or false to decide whether to run Pattern search (Pattern search) 276stg.psearch = false; 277 278% Options for Pattern search (Pattern search options) 279stg.psearch_options = optimoptions(@patternsearch,... 280 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 281 'UseCompletePoll',true,'UseCompleteSearch',true,... 282 'MaxMeshSize',2,'MaxFunctionEvaluations',2000); 283 284% True or false to decide whether to run Genetic algorithm (Genetic 285% algorithm) 286stg.ga = true; 287 288% Options for Genetic algorithm (Genetic algorithm options) 289stg.ga_options = optimoptions(@ga,'MaxGenerations',200,... 290 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 291 'PopulationSize',stg.popsize,... 292 'MutationFcn','mutationadaptfeasible','Display','diagnose'); 293 294% True or false to decide whether to run Particle swarm (Particle swarm) 295stg.pswarm = false; 296 297% Options for Particle swarm (Particle swarm options) 298stg.pswarm_options = optimoptions('particleswarm',... 299 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 300 'SwarmSize',stg.popsize); 301 302% True or false to decide whether to run Surrogate optimization (Surrogate 303% optimization) 304stg.sopt = false; 305 306% Options for Surrogate optimization (Surrogate optimization options) 307stg.sopt_options = optimoptions('surrogateopt',... 308 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 309 'MaxFunctionEvaluations',5000,... 310 'MinSampleDistance',0.2,'MinSurrogatePoints',32*2+1); 311end
Import
Default settings code
1% True or false to decide whether to run import functions (Import) 2stg.import = true; 3 4 5% Name of the excel file with the sbtab (SBtab excel name) 6stg.sbtab_excel_name = "SBTAB.xlsx"; 7 8% Name of the model (Name) 9stg.name = "model_name"; 10 11% Name of the default model compartment (Compartment name) 12stg.cname = "Compartment"; 13 14% Name of the sbtab saved in .mat format (SBtab name) 15stg.sbtab_name = "SBtab_" + stg.name;Example settings code
1% True or false to decide whether to run import functions (Import) 2stg.import = true; 3 4% Name of the excel file with the sbtab (SBtab excel name) 5stg.sbtab_excel_name = "SBTAB example.xlsx"; 6 7% Name of the model (Name) 8stg.name = "Example"; 9 10% Name of the default model compartment (Compartment name) 11stg.cname = "Compartment";
stg.import - (logical) Decide whether to run import functions
stg.sbtab_excel_name - (string) Name of the Excel file with the SBtab
stg.name - (string) Name of the model
stg.cname - (string) Name of the default model compartment
stg.sbtab_name - (string) Name of the SBtab saved in .mat format
Analysis
Default settings code
1% Experiments to run (example experiment 1 to 3 and experimet 6) 2stg.exprun = [1:3,6]; 3 4% Choice between 0,1,2 and 3,4 to choose the scoring method (check 5% documentation) (Use logarithm) 6stg.useLog = 4; 7 8% True or false to decide whether to use multicore everywhere it is 9% available (Optimization Multicore) 10stg.optmc = true; 11 12% Choice of ramdom seed (Ramdom seed) 13stg.rseed = 1; 14 15% True or false to decide whether to use display simulation diagnostics in 16% the console (Simulation Console) 17stg.simcsl = false; 18 19% True or false to decide whether to display optimization results on 20% console (Optimization console) 21stg.optcsl = false; 22 23% True or false to decide whether to save results (Save results) 24stg.save_results = true; 25 26% True or false to decide whether to run detailed simulation for plots 27stg.simdetail = true;Example settings code
1% Experiments to run 2stg.exprun = [1,2]; 3% stg.exprun = [1,2,3]; 4 5% Choice between 0,1,2 and 3 to change either and how to apply log10 to the 6% scores (check documentation) (Use logarithm) 7stg.useLog = 0; 8 9% True or false to decide whether to use multicore everywhere it is 10% available (Optimization Multicore) 11stg.optmc = false; 12 13% Choice of ramdom seed (Ramdom seed) 14stg.rseed = 1; 15 16% True or false to decide whether to use display simulation diagnostics in 17% the console (Simulation Console) 18stg.simcsl = false; 19 20% True or false to decide whether to display optimization results on 21% console (Optimization console) 22stg.optcsl = true; 23 24% True or false to decide whether to display PLA results on console (PLA 25% console) 26stg.placsl = true; 27 28% True or false to decide whether to save results (Save results) 29stg.save_results = true; 30 31% True or false to decide whether to run detailed simulation for plots 32stg.simdetail = false;
stg.exprun - (double) Experiments to run
stg.useLog - (double) Choice between 0,1,2 and 3 to change either and how to apply log10 to the scores, check results:
stg.optmc - (logical) Decide whether to use multicore everywhere it is available
stg.rseed - (double) Choice of random seed
stg.simcsl - (logical) Decide whether to display simulation diagnostics in the console
stg.optcsl - (logical) Decide whether to display optimization results on console
stg.save_results - (logical) Decide whether to save results
stg.simdetail - (logical) Decide whether to run detailed simulation for plots
Simulation
Default settings code
1% Maximum time for each individual function to run in seconds (Maximum 2% time) 3stg.maxt = 10; 4 5% Equilibration time (Equilibration time) 6stg.eqt = 50000; 7 8% True or false to decide whether to do Dimensional Analysis (Dimensional 9% Analysis) 10stg.dimenanal = false; 11 12% True or false to decide whether to do Unit conversion (Unit conversion) 13stg.UnitConversion = false; 14 15% True or false to decide whether to do Absolute Tolerance Scaling 16% (Absolute Tolerance Scaling) 17stg.abstolscale = true; 18 19% Value of Relative tolerance (Relative tolerance) 20stg.reltol = 1.0E-4; 21 22% Value of Absolute tolerance (Absolute tolerance) 23stg.abstol = 1.0E-4; 24 25% Time units for simulation (Simulation time) 26stg.simtime = "second"; 27 28% True or false to decide whether to run sbioaccelerate (after changing 29% this value you need to run "clear functions" to see an effect) 30% (sbioaccelerate) 31stg.sbioacc = true; 32 33% (Absolute tolerance step size for equilibration) 34stg.abstolstepsize_eq = []; 35 36% Max step size in the simulation (if empty matlab decides whats best) 37% (Maximum step) 38stg.maxstep = [10]; 39 40% Max step size in the equilibration (if empty matlab decides whats best) 41% (Maximum step) 42stg.maxstepeq = [2]; 43 44% Max step size in the detailed plots (if empty matlab decides whats best) 45% (Maximum step) 46stg.maxstepdetail = [2]; 47 48% Default score when there is a simulation error, this is needed to keep 49% the optimizations working. 50% (error score) 51stg.errorscore = 10^5;Example settings code
1% Maximum time for each individual function to run in seconds (Maximum 2% time) 3stg.maxt = 2; 4 5% Equilibration time (Equilibration time) 6stg.eqt = 50000; 7 8% True or false to decide whether to do Dimensional Analysis (Dimensional 9% Analysis) 10stg.dimenanal = true; 11 12% True or false to decide whether to do Unit conversion (Unit conversion) 13stg.UnitConversion = true; 14 15% True or false to decide whether to do Absolute Tolerance Scaling 16% (Absolute Tolerance Scaling) 17stg.abstolscale = true; 18 19% Value of Relative tolerance (Relative tolerance) 20stg.reltol = 1.0E-4; 21 22% Value of Absolute tolerance (Absolute tolerance) 23stg.abstol = 1.0E-7; 24 25% Time units for simulation (Simulation time) 26stg.simtime = "second"; 27 28% True or false to decide whether to run sbioaccelerate (after changing 29% this value you need to run "clear functions" to see an effect) 30% (sbioaccelerate) 31stg.sbioacc = false; 32 33% Max step size in the simulation (if empty matlab decides whats best) 34% (Maximum step) 35stg.maxstep = []; 36 37% Max step size in the equilibration (if empty matlab decides whats best) 38% (Maximum step) 39stg.maxstepeq = []; 40 41% Max step size in the detailed plots (if empty matlab decides whats best) 42% (Maximum step) 43stg.maxstepdetail = [0.001]; 44 45% Default score when there is a simulation error, this is needed to keep 46% the optimizations working. (error score) 47stg.errorscore = 10^5;
stg.maxt - (double) Maximum time for each individual function has to run in seconds
stg.eqt - (double) Equilibration time in seconds
stg.dimenanal - (logical) Decide whether to do Dimensional Analysis
stg.UnitConversion - (logical) Decide whether to do Unit conversion
stg.abstolscale - (logical) Decide whether to do Absolute Tolerance Scaling
stg.reltol - (double) Value of Relative tolerance
stg.abstol - (double) Value of Absolute tolerance
stg.simtime - (string) Time units for simulation
stg.sbioacc - (logical) Decide whether to run sbioaccelerate (after changing this value you need to run “clear functions” to see an effect)
stg.abstolstepsize_eq - (double) Absolute tolerance step size for equilibration (if empty MATLAB® decides whats best)
stg.maxstep - (double) Max step size in the simulation (if empty MATLAB® decides what’s best)
stg.maxstepeq - (double) Max step size in the equilibration (if empty MATLAB® decides whats best)
stg.maxstepdetail - (double) Max step size in the detailed plots (if empty MATLAB® decides whats best)
stg.errorscore - (double) Default score when there is a simulation error, this is needed to keep the optimizations working.
Model
Default settings code
1% Number of parameters to optimize (Parameter number) 2stg.parnum = 5; 3 4original_parameter_set = zeros(1,10); 5 6% Array with the lower bound of all parameters (Lower bound) 7stg.lb = original_parameter_set-5; 8 9% Array with the upper bound of all parameters (Upper bound) 10stg.ub = original_parameter_set+5;Example settings code
1% Number of parameters to optimize (Parameter number) 2stg.parnum = 12; 3 4% Index for the parameters that have thermodynamic constrains (Termodiamic 5% Constrains Index) 6stg.tci = [8]; 7 8% Parameters to multiply to the first parameter (in Stg.partest to get to 9% the correct thermodynamic constrain formula) (Termodiamic Constrains 10% multipliers) 11stg.tcm([8],1) = [4]; 12stg.tcm([8],2) = [5]; 13stg.tcm([8],3) = [7]; 14 15% Parameters to divide to the first parameter (in Stg.partest to get to the 16% correct thermodynamic constrain formula) (Termodiamic Constrains 17% divisors) 18stg.tcd([8],1) = [1]; 19stg.tcd([8],2) = [3]; 20stg.tcd([8],3) = [6]; 21 22% Array with the lower bound of all parameters (Lower bound) 23stg.lb = zeros(1,stg.parnum)-4; 24 25% Array with the upper bound of all parameters (Upper bound) 26stg.ub = zeros(1,stg.parnum)+4;
stg.parnum - (double) Number of parameters to optimize
stg.tci - (double) Index for the parameters that have thermodynamic constraints
stg.tcm - (double) Parameters to multiply to the first parameter (in stg.partest to get to the correct thermodynamic constraint formula)
stg.tcd - (double) Parameters to divide to the first parameter (in stg.partest to get to the correct thermodynamic constraint formula)
stg.lb - (double) Lower bound of all parameters
\[stg.lb = \begin{bmatrix} lb_{1} & lb_{2} & ... & lb_{i} \end{bmatrix}\]\(i =\) Parameter index
stg.ub - (double) Upper bound of all parameters
\[stg.up = \begin{bmatrix} ub_{1} & ub_{2} & ... & ub_{i} \end{bmatrix}\]\(i =\) Parameter index
Diagnostics
Default settings code
1% Choice of what parameters in the array to test, the indices correspond to 2% the parameters in the model and the numbers correspond to the parameters 3% in the optimization array, usually not all parameters are optimized so 4% there needs to be a match between one and the other. (Parameters to test) 5% In this example there are ten parameters in this imaginary model and we 6% are only interested in parameter 2,4,8,9, and 10. Note that stg.parnum is 7% five because of this and not ten 8stg.partest(:,1) = [0,1,0,2,0,0,0,3,4,5]; 9 10% (Parameter array to test) 11stg.pat = [1:2]; 12 13% All the parameter arrays, in this case there is only one (parameters here 14% are in log10 space)(Parameter arrays) 15stg.pa(1,:) = [1,1,1,1,1]; 16stg.pa(1,:) = [1,0,1,2,1]; 17 18% Best parameter array found so far for the model (Best parameter array) 19stg.bestpa = stg.pa(1,:);Example settings code
1% Choice of what parameters in the array to test, the indices correspond to 2% the parameters in the model and the numbers correspond to the parameters 3% in the optimization array, usually not all parameters are optimized so 4% there needs to be a match between one and the other. (Parameters to test) 5stg.partest(:,1) = [1 ,2 ,3 ,4 ,5 ,6 ,7 ,2 ,8 ,9,... 6 10 ,11 ,12]; 7 8% (Parameter array to test) 9stg.pat = 1:3; 10 11% All the parameter arrays, in this case there is only one (Parameter 12% arrays) 13stg.pa(1,:) = [3.999,0.696,1.072,3.429,-0.751,-3.741,-0.569,0.831,3.068,0.921,-2.156,-1.970]; 14stg.pa(2,:) = stg.pa(1,:)-1; 15stg.pa(3,:) = stg.pa(1,:)+1; 16 17% Best parameter array found so far for the model (Best parameter array) 18stg.bestpa = stg.pa(1,:);
stg.partest - (double) Choice of which parameters to work on, since depending on the task, not all SBtab parameters are worked on. k indices correspond to the parameters in the SBtab and numbers up to i correspond to the parameters in the work set. This is the set that actually gets used for diagnostics, optimization, and sensitivity analyis.
\[stg.partest_k = \begin{bmatrix} 1_{k_1} & 2_{k_2} & ... & i_{k_{end}} \end{bmatrix}\]In our example model parameter 216 from the SBtab is parameter number 1 of the work set, parameter 217 from the SBtab is parameter number 2 of the work set, and successively.
\[stg.partest_{[216:227]} = \begin{bmatrix} 1_{216} & 2_{217} & ... & 6_{221} & 1_{222} & 2_{223} & ... & 6_{227} \end{bmatrix}\]stg.pat - (double) Index(\(j\)) of the parameter set to work on
stg.pa - (double) All the parameter sets
\[\begin{split}stg.pa = \begin{bmatrix} x_{1,1} & x_{2,1} & ... & x_{i,1} \\ x_{1,2} & x_{2,2} & ... & x_{i,2} \\ ... & ... & ... & ... \\ x_{1,j} & x_{2,j} & ... & x_{i,j} \end{bmatrix}\end{split}\]stg.bestpa - (double) Best parameter set found so far during optimization
\[stg.bestx = \begin{bmatrix} bestx_{1} & bestx_{2} & ... & bestx_{i} \end{bmatrix}\]\(x =\) Parameters being worked on
\(i =\) Index of Parameters being worked on
\(k =\) Index of the parameters in SBtab
\(j =\) Index of the Parameter set to work on
Plots
Default settings code
1% True or false to decide whether to plot results (Plots) 2stg.plot = true; 3 4% True or false to decide whether to use long names in the title of the 5% outputs plots in f_plot_outputs.m (Plot outputs long names) 6stg.plotoln = true;Example settings code
1% True or false to decide whether to plot results (Plots) 2stg.plot = true; 3 4% True or false to decide whether to use long names in the title of the 5% outputs plots in f_plot_outputs.m (Plot outputs long names) 6stg.plotnames = true;
stg.plot - (logical) Decide whether to plot results
stg.plotoln - (logical) Decide whether to use long names in the title of the output plots in f_plot_outputs.m
Global Sensitivity Analysis (GSA)
Default settings code
1% Number of samples to use in SA (Sensitivity analysis number of samples) 2stg.sansamples = 100; 3 4% True or false to decide whether to subtract the mean before calculating 5% SI and SIT (Sensitivity analysis subtract mean) 6stg.sasubmean = true; 7 8% Choose the way you want to obtain the samples of the parameters for 9% performing the SA; 0 Log uniform distribution truncated at the parameter 10% bounds 1 Log normal distribution with mu as the best value for a 11% parameter and sigma as stg.sasamplesigma truncated at the parameter 12% bounds 2 same as 1 without truncation 3 Log normal distribution centered 13% at the mean of the parameter bounds and sigma as stg.sasamplesigma 14% truncated at the parameter bounds 4 same as 3 without truncation. 15% (Sensitivity analysis sampling mode) 16stg.sasamplemode = 2; 17 18% Sigma for creating the normal distribution of parameters to perform 19% sensitivity analysis (Sensitivity analysis sampling sigma) 20stg.sasamplesigma = 0.1;Example settings code
1% Number of samples to use in SA (Sensitivity analysis number of samples) 2stg.sansamples = 1000; 3 4% True or false to decide whether to subtract the mean before calculating 5% SI and SIT (Sensitivity analysis subtract mean) 6stg.sasubmean = true; 7 8% Choose the way you want to obtain the samples of the parameters for 9% performing the SA; 0 Log uniform distribution truncated at the parameter 10% bounds 1 Log normal distribution with mu as the best value for a 11% parameter and sigma as stg.sasamplesigma truncated at the parameter 12% bounds 2 same as 1 without truncation 3 Log normal distribution centered 13% at the mean of the parameter bounds and sigma as stg.sasamplesigma 14% truncated at the parameter bounds 4 same as 3 without truncation. 15% (Sensitivity analysis sampling mode) 16stg.sasamplemode = 2; 17 18% Sigma for creating the normal distribution of parameters to perform 19% sensitivity analysis (Sensitivity analysis sampling sigma) 20stg.sasamplesigma = 0.1; 21 22stg.gsabootstrapsize = ceil(sqrt(stg.sansamples)); 23 24%% Profile Likelihood 25 26% Parameter(optimization array) that is being worked on in a specific 27% iteration of PL (if -1 no parameter is being worked in PL) 28% (Profile Likelihood Index) 29stg.PLind = -1; 30 31% Which parameters to do PL on, it should be all parameters but can also be 32% a subset for testing purposes 33% (Profile Likelihood parameters to Test) 34stg.pltest = (1:12); 35 36% How many points to do for each parameter in the PL 37% (Profile Likelihood Resolution) 38stg.plres = 80; 39 40% True or false to decide whether to do plots after calculating PL 41% (Profile Likelihood Plots) 42stg.plplot = true; 43 44% True or false to decide whether to run simulated annealing 45% (Profile Likelihood Simulated Annealing) 46stg.plsa = true; 47 48% Options for simulated annealing 49stg.plsao = optimoptions(@simulannealbnd,'Display','off', ... 50 'InitialTemperature',... 51 ones(1,stg.parnum)*1,'MaxTime',5,'ReannealInterval',40); 52 53% 0 or 1 to decide whether to run fmincon 54% (Profile Likelihood FMincon) 55stg.plfm = false; 56 57% Options for fmincon 58stg.plfmo = optimoptions('fmincon','Display','off',... 59 'Algorithm','interior-point',... 60 'MaxIterations',5,'OptimalityTolerance',0,... 61 'StepTolerance',1e-6,'FiniteDifferenceType','central');
stg.sansamples - (double) Number of samples to use in GSA (in total (2+npars)*sansamples simulations will be performed, where npars are the number of parameters).
stg.sasubmean - (logical) Decide whether to subtract mean before calculating SI and STI, see Halnes et al 2009.
stg.sasamplemode - (double) Choose the way you want to obtain the samples of the parameters for performing the GSA;
Reciprocal (log uniform) distribution
\(X_{i} \sim Reciprocal(a_{i},b_{i})\)
\(i =\) Parameter index
\(a_{i} = stg.lb_{i}\)
\(b_{i} = stg.ub_{i}\)
Example distribution with \(a = -1, b = 1\)
Log normal distribution with μ corresponding to the best value for a parameter, as recieved from the optimization, and σ as stg.sasamplesigma truncated at the parameter bounds
\(X_{i} \sim TruncatedLogNormal(μ_{i}, σ, a_{i}, b_{i})\)
\(i =\) Parameter index
\(μ_{i} = bestx_{i}\)
\(σ = stg.sasamplesigma\)
\(a_{i} = stg.lb_{i}\)
\(b_{i} = stg.ub_{i}\)
Example distribution with \(μ = 0.5, σ = 1, a = -1, b = 1\)
same as 1 without truncation
\(X_{i} \sim LogNormal(μ, σ)\)
\(i =\) Parameter index
\(μ_{i} = bestx_{i}\)
\(σ = stg.sasamplesigma\)
Example distribution with \(μ = 0.5, σ = 1\)
Log normal distribution with μ corresponding to the mean of the parameter bounds and σ as stg.sasamplesigma but truncated at the parameter bounds
\(X_{i} \sim TruncatedLogNormal(μ_{i}, σ, a_{i}, b_{i})\)
\(i =\) Parameter index
\(μ_{i} = \frac{stg.lb_{i} + (stg.ub_{i} - stg.lb_{i})}{2}\)
\(σ = stg.sasamplesigma\)
\(a_{i} = stg.lb_{i}\)
\(b_{i} = stg.ub_{i}\)
Example distribution with \(μ = \frac{a+(b-a)}{2}, σ = 1, a = -1, b = 1\)
same as 3 without truncation.
\(X_{i} \sim LogNormal(mu_{i}, σ)\)
\(i =\) Parameter index
\(μ_{i} = \frac{stg.lb_{i} + (stg.ub_{i} - stg.lb_{i})}{2}\)
\(σ = stg.sasamplesigma\)
Example distribution with \(μ = \frac{a+(b-a)}{2}, σ = 1, a = -1, b = 1\)
stg.sasamplesigma - (double) σ for creating the normal distribution of parameters to perform sensitivity analysis
Optimization
Default settings code
1% Time for the optimization in seconds (fmincon does not respect this 2% time!!) (Optimization time) 3stg.optt = 60*5; 4 5% Population size for the algorithms that use populations (Population size) 6stg.popsize = 144; 7 8% optimization start method, choose between: 1 Random starting point or 9% group of starting points inside the bounds 2 Random starting point or 10% group of starting points near the best point (Optimization start method) 11stg.osm = 1; 12 13% Distance from best parameter array to be used in stg.osm method 2 14% (Distance from best parameter array) 15stg.dbs = 0.1; 16 17% True or false to decide whether to use Multistart (Multistart) 18stg.mst = false; 19 20% Multistart size 21stg.msts = 1; 22 23% True or false to decide whether to display Plots (Plots doesn't work if 24% using multicore) (Optimization plots) 25stg.optplots = true; 26 27% True or false to decide whether to run fmincon (no gradient so this 28% doesn't work very well, no max time!!) 29stg.fmincon = false; 30 31% Options for fmincon (fmincon options) 32stg.fm_options = optimoptions('fmincon',... 33 'Algorithm','interior-point',... 34 'MaxIterations',2,'OptimalityTolerance',0,... 35 'StepTolerance',1e-6,'FiniteDifferenceType','central',... 36 'MaxFunctionEvaluations',10000); 37 38% True or false to decide whether to run simulated annealing (Simulated 39% annealing) 40stg.sa = false; 41 42% Options for simulated annealing (Simulated annealing options) 43stg.sa_options = optimoptions(@simulannealbnd, ... 44 'MaxTime',stg.optt,... 45 'ReannealInterval',40); 46 47% True or false to decide whether to run Pattern search (Pattern search) 48stg.psearch = false; 49 50% Options for Pattern search (Pattern search options) 51%stg.psearch_options = optimoptions(@patternsearch,... 52 %'MaxTime',stg.optt,'UseParallel',stg.optmc,... 53 %'UseCompletePoll',true,'UseCompleteSearch',true,... 54 %'MaxMeshSize',2,'MaxFunctionEvaluations',2000); 55 56% True or false to decide whether to run Genetic algorithm (Genetic 57% algorithm) 58stg.ga = false; 59 60% Options for Genetic algorithm (Genetic algorithm options) 61stg.ga_options = optimoptions(@ga,'MaxGenerations',200,... 62 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 63 'PopulationSize',stg.popsize,... 64 'MutationFcn','mutationadaptfeasible'); 65 66% True or false to decide whether to run Particle swarm (Particle swarm) 67stg.pswarm = false; 68 69% Options for Particle swarm (Particle swarm options) 70stg.pswarm_options = optimoptions('particleswarm',... 71 'MaxTime',stg.optt,'UseParallel',stg.optmc,'MaxIterations',200,... 72 'SwarmSize',stg.popsize); 73 74% True or false to decide whether to run Surrogate optimization (Surrogate 75% optimization) 76stg.sopt = false; 77 78% Options for Surrogate optimization (Surrogate optimization options) 79stg.sopt_options = optimoptions('surrogateopt',... 80 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 81 'MaxFunctionEvaluations',5000,... 82 'MinSampleDistance',0.2,'MinSurrogatePoints',32*2+1); 83endExample settings code
1% Time for the optimization in seconds (fmincon does not respect this 2% time!!) (Optimization time) 3stg.optt = 60*1; 4 5% Population size for the algorithms that use populations (Population size) 6stg.popsize = 10; 7 8% optimization start method, choose between: 1 Random starting point or 9% group of starting points inside the bounds 2 Random starting point or 10% group of starting points near the best point (Optimization start method) 11stg.osm = 1; 12 13% Distance from best parameter array to be used in stg.osm method 2 14% (Distance from best parameter array) 15stg.dbs = 0.1; 16 17% True or false to decide whether to use Multistart (Multistart) 18stg.mst = false; 19 20% Multistart size 21stg.msts = 1; 22 23% True or false to decide whether to display Plots (Plots doesn't work if 24% using multicore) (Optimization plots) 25stg.optplots = true; 26 27% True or false to decide whether to run fmincon (no gradient so this 28% doesn't work very well, no max time!!) 29stg.fmincon = false; 30 31% Options for fmincon (fmincon options) 32stg.fm_options = optimoptions('fmincon',... 33 'UseParallel',stg.optmc,... 34 'Algorithm','interior-point',... 35 'MaxIterations',2,'OptimalityTolerance',0,... 36 'StepTolerance',1e-6,'FiniteDifferenceType','central',... 37 'MaxFunctionEvaluations',10000); 38 39% True or false to decide whether to run simulated annealing (Simulated 40% annealing) 41stg.sa = false; 42 43% Options for simulated annealing (Simulated annealing options) 44stg.sa_options = optimoptions(@simulannealbnd, ... 45 'MaxTime',stg.optt,... 46 'InitialTemperature',... 47 ones(1,stg.parnum)*2,'ReannealInterval',40); 48 49% True or false to decide whether to run Pattern search (Pattern search) 50stg.psearch = false; 51 52% Options for Pattern search (Pattern search options) 53stg.psearch_options = optimoptions(@patternsearch,... 54 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 55 'UseCompletePoll',true,'UseCompleteSearch',true,... 56 'MaxMeshSize',2,'MaxFunctionEvaluations',2000); 57 58% True or false to decide whether to run Genetic algorithm (Genetic 59% algorithm) 60stg.ga = true; 61 62% Options for Genetic algorithm (Genetic algorithm options) 63stg.ga_options = optimoptions(@ga,'MaxGenerations',200,... 64 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 65 'PopulationSize',stg.popsize,... 66 'MutationFcn','mutationadaptfeasible','Display','diagnose'); 67 68% True or false to decide whether to run Particle swarm (Particle swarm) 69stg.pswarm = false; 70 71% Options for Particle swarm (Particle swarm options) 72stg.pswarm_options = optimoptions('particleswarm',... 73 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 74 'SwarmSize',stg.popsize); 75 76% True or false to decide whether to run Surrogate optimization (Surrogate 77% optimization) 78stg.sopt = false; 79 80% Options for Surrogate optimization (Surrogate optimization options) 81stg.sopt_options = optimoptions('surrogateopt',... 82 'MaxTime',stg.optt,'UseParallel',stg.optmc,... 83 'MaxFunctionEvaluations',5000,... 84 'MinSampleDistance',0.2,'MinSurrogatePoints',32*2+1); 85end
stg.optt - (double) Time for the optimization in seconds (fmincon does not respect this time!!)
stg.popsize - (double) Population size (for the algorithms that use populations)
stg.osm - (double) optimization start method, choose between
Get a random starting parameter set or group of starting parameter sets inside the bounds
Get a random starting parameter set or group of starting parameter sets near the best parameter set
stg.dbpa - (double) Distance from best parameter set to be used in stg.osm method 2
stg.mst - (logical) Decide whether to use one or multiple starting parameter sets for the optimization
stg.msts - (double) Number of starting parameter sets for the optimizations
stg.optplots - (logical) Decide whether to display optimiazation plots (They aren’t ploted if running the code in multicore)
stg.fmincon - (logical) Decide whether to run fmincon (no gradient in our models so this doesn’t work very well, does not respect time set for the optimization!!)
stg.fm_options - (optim.options.Fmincon) Options for fmincon
stg.sa - (logical) Decide whether to run simulated annealing
stg.sa_options - (optim.options.SimulannealbndOptions) Options for simulated annealing
stg.psearch - (logical) Decide whether to run Pattern search
stg.psearch_options - (optim.options.PatternsearchOptions) Options for Pattern search
stg.ga - (logical) Decide whether to run Genetic algorithm
stg.ga_options - (optim.options.GaOptions) Options for Genetic algorithm
stg.pswarm - (logical) Decide whether to run Particle swarm
stg.pswarm_options - (optim.options.Particleswarm) Options for Particle swarm
stg.sopt - (logical) Decide whether to run Surrogate optimization
stg.sopt_options - (optim.options.Surrogateopt) Options for Surrogate optimization
Automatically generated at Import
stg.expn - (double) Total number of experiments stored in the SBtab
stg.outn - (double) Total number of experimental outputs specified in the SBtab
References
Results
Scoring and saved simulation output
Every time a simulation is run the simulated results are compared to the results provided and a score is calculated. Additionally the end point of the experimental output of all simulations is also stored. When performing the diagnostics function an MATLAB® representation of the entire run is also saved.
simd - Simulation results (MATLAB® representation)
st - Total score
To simplify representations the following correspondence has been used
\(score_{i,j,k} = \frac{1}{n} \sum_{i=1}^n \left(\frac{Y_{i,j,k}-y_{i,j,k}}{τ_{i,j,k}}\right)^2\)
if stg.useLog = 0
\(st(θ;Y,τ) = \sum_{k=1}^l \sum_{j=1}^m score_{i,j,k}\)
if stg.useLog = 1
\(st(θ;Y,τ) = \sum_{k=1}^l \sum_{j=1}^m log_{10}(score_{i,j,k})\)
if stg.useLog = 2
\(st(θ;Y,τ) = \sum_{k=1}^l log_{10}(\sum_{j=1}^m score_{i,j,k})\)
if stg.useLog = 3
\(st(θ;Y,τ) = log_{10}(\sum_{k=1}^l \sum_{j=1}^m score_{i,j,k})\)
se - Score of each experiment
if stg.useLog = 0 or 3
\[\begin{split}se(θ;Y,τ) = \begin{bmatrix} \sum_{j=1}^m score_{i,j,1} \\ \sum_{j=1}^m score_{i,j,2} \\ ... \\ \sum_{j=1}^m score_{i,j,k} \end{bmatrix}\end{split}\]if stg.useLog = 1
\[\begin{split}se(θ;Y,τ) = \begin{bmatrix} \sum_{j=1}^m log_{10}(score_{i,j,1}) \\ \sum_{j=1}^m log_{10}(score_{i,j,2}) \\ ... \\ \sum_{j=1}^m log_{10}(score_{i,j,k}) \end{bmatrix}\end{split}\]if stg.useLog = 2
\[\begin{split}se(θ;Y,τ) = \begin{bmatrix} log_{10}(\sum_{j=1}^m score_{i,j,1}) \\ log_{10}(\sum_{j=1}^m score_{i,j,2})\\ ... \\ log_{10}(\sum_{j=1}^m score_{i,j,k}) \end{bmatrix}\end{split}\]sd - Score of each experimental outputs in all experiments
if stg.useLog = 0,2, or 3
\[\begin{split}sd(θ;Y,τ) = \begin{bmatrix} score_{i,1,1} & score_{i,2,1} & ... & score_{i,j,1}\\ score_{i,1,2} & score_{i,2,2} & ... & score_{i,j,2}\\ ... & ... & ... & ... \\ score_{i,1,k} & score_{i,2,k} & ... & score_{i,j,k} \end{bmatrix}\end{split}\]if stg.useLog = 1
\[\begin{split}sd(θ;Y,τ) = \begin{bmatrix} log_{10}(score_{i,1,1}) & log_{10}(score_{i,2,1}) & ... & log_{10}(score_{i,j,1})\\ log_{10}(score_{i,1,2}) & log_{10}(score_{i,2,2}) & ... & log_{10}(score_{i,j,2})\\ ... & ... & ... & ... \\ log_{10}(score_{i,1,k}) & log_{10}(score_{i,2,k}) & ... & log_{10}(score_{i,j,k}) \end{bmatrix}\end{split}\]xfinal - Value of each experimental outputs at the end of the simulation
\[\begin{split}xfinal(θ;Y,τ) = \begin{bmatrix} y_{n,1,1} & y_{n,2,1} & ... & y_{n,j,1} \\ y_{n,1,2} & y_{n,2,2} & ... & y_{n,j,2} \\ ... & ... & ... & ... \\ y_{n,1,k} & y_{n,2,k} & ... & y_{n,j,k} \end{bmatrix}\end{split}\]\(F =\) Objective function for Particle Swarm optimization
\(Y =\) Data provided for fitting
\(y =\) Simulation results of the updated model under parameterization \(θ\)
\(θ =\) New parameterization for \(y\)
\(τ =\) Allowed mismatch between the two simulation results, analogous to the standard deviation of a Gaussian noise model in data fitting
\(n/i =\) Number/index of points in a given experimental output
\(m/j =\) Number/index of experimental outputs
\(l/k =\) Number/index of experiments
Diagnostics
When running the diagnostics a struct gets created that stores all the oputputs of the f_sim_score function.
rst.diag.simd - Simulation results (MATLAB® representation)
rst.diag.st - Total score
rst.diag.se - Score per experiment
rst.diag.sd - Score per experimental outputs in all experiments
rst.diag.xfinal - x value of all the species being tested at the end of the simulation
Optimization
rst.opt.name - Name of optimizer that was used
rst.opt.x - Best parameter set found by the optimization
rst.opt.fval - Score for that best parameter set
rst.opt.exitflag - Diagnostics to see how the optimization went
rst.opt.output - Diagnostics to see how the optimization went
Sensitivity Analysis
The calculations performed to obtain these sensitivities where performed according to the equations described in Halnes et al 2009.
rst.SA.M1 - Matrix with (\(r*k\)) random numbers within the lower and upper bound ranges set for each parameter
\[\begin{split}M_1 = \begin{bmatrix} x_{1}^{(1)} & x_{2}^{(1)} & ... & x_{k}^{(1)} \\ x_{1}^{(2)} & x_{2}^{(2)} & ... & x_{k}^{(2)} \\ ... & ... & ... & ... \\ x_{1}^{(r)} & x_{2}^{(r)} & ... & x_{k}^{(r)} \end{bmatrix}\end{split}\]\(x =\) Parameters
\(k =\) Total number of parameters (stg.parnum)
\(r =\) Total number of Samples (stg.sansamples)
rst.SA.M2 - Same as rst.SA.M1 but different random initialization
\[\begin{split}M_2 = \begin{bmatrix} x_{1}^{(1')} & x_{2}^{(1')} & ... & x_{k}^{(1')} \\ x_{1}^{(2')} & x_{2}^{(2')} & ... & x_{k}^{(2')} \\ ... & ... & ... & ... \\ x_{1}^{(r')} & x_{2}^{(r')} & ... & x_{k}^{(r')} \end{bmatrix}\end{split}\]\(x =\) Parameters
\(k =\) Total number of parameters (stg.parnum)
\(r =\) Total number of Samples (stg.sansamples)
rst.SA.N - Matrix of size (\(r*k*k\)) with columns exchanged between M1 and M2 as follows:
\[\begin{split}N_i = \begin{bmatrix} x_{1}^{(1')} & x_{2}^{(1')} & ... & x_{i}^{(1)} & ... & x_{k}^{(1')} \\ x_{1}^{(2')} & x_{2}^{(2')} & ... & x_{i}^{(2)} & ... & x_{k}^{(2')} \\ ... & ... & ... & ... & ... & ... \\ x_{1}^{(r')} & x_{2}^{(r')} & ... & x_{i}^{(r)} & ... & x_{k}^{(r')} \end{bmatrix}\end{split}\]\(x =\) Parameters
\(k =\) Total number of parameters (stg.parnum)
\(r =\) Total number of Samples (stg.sansamples)
\(i =\) Index of each parameter
rst.SA.fM1 -
\[\begin{split}fM_1 = \begin{bmatrix} f(M_1^{(1)}) \\ f(M_1^{(2)}) \\ ... \\ f(M_1^{(r)}) \end{bmatrix} = \begin{bmatrix} f(x_{1}^{(1)} & x_{2}^{(1)} & ... & x_{k}^{(1)}) \\ f(x_{1}^{(2)} & x_{2}^{(2)} & ... & x_{k}^{(2)}) \\ ... & ... & ... & ... \\ f(x_{1}^{(r)} & x_{2}^{(r)} & ... & x_{k}^{(r)}) \end{bmatrix}\end{split}\]\(k =\) Total number of parameters (stg.parnum)
\(r =\) Total number of Samples (stg.sansamples)
rst.SA.fM2 -
\[\begin{split}fM_2 = \begin{bmatrix} f(M_2^{(1')}) \\ f(M_2^{(2')}) \\ ... \\ f(M_2^{(r')}) \end{bmatrix} = \begin{bmatrix} f(x_{1}^{(1')} & x_{2}^{(1')} & ... & x_{k}^{(1')}) \\ f(x_{1}^{(2')} & x_{2}^{(2')} & ... & x_{k}^{(2')}) \\ ... & ... & ... & ... \\ f(x_{1}^{(r')} & x_{2}^{(r')} & ... & x_{k}^{(r')}) \end{bmatrix}\end{split}\]\(k =\) Total number of parameters (stg.parnum)
\(r =\) Total number of Samples (stg.sansamples)
rst.SA.fN -
\[\begin{split}fN_i = \begin{bmatrix} f(N_i^{(1)}) \\ f(N_i^{(2)}) \\ ... \\ f(N_i^{(r)}) \end{bmatrix} = \begin{bmatrix} f(x_{1}^{(1')} & x_{2}^{(1')} & ... & x_{i}^{(1)} & ... & x_{k}^{(1')}) \\ f(x_{1}^{(2')} & x_{2}^{(2')} & ... & x_{i}^{(2)} & ... & x_{k}^{(2')}) \\ ... & ... & ... & ... & ... & ... \\ f(x_{1}^{(r')} & x_{2}^{(r')} & ... & x_{i}^{(r)} & ... & x_{k}^{(r')}) \end{bmatrix}\end{split}\]\(k =\) Total number of parameters (stg.parnum)
\(r =\) Total number of Samples (stg.sansamples)
\(i =\) Index of each parameter
rst.SA.SI - First order effects
\(S_{i}=\frac{V_{Θ_{i}}(E_{Θ_{-i}}(Y|Θ_{i}))}{V(Y)}=\frac{U_{i}-E^2(Y)}{V(Y)}\)
\(U_{i}=\frac{1}{n-1}\sum_{r=1}^nf(M_1^r)f(N_i^r)\)
\(E^2(Y)=\frac{1}{n}\sum_{r=1}^nf(M_1^r)f(M_2^r)\)
\(V(Y) = \frac{1}{n-1}f^2(M_1^r)-E^2(Y)\)
\(V =\) Variance
\(E(... |...) =\) Conditional expected value
\(Θ =\) Parameters of the model
\(Y =\) Scalar output from the model
\(n =\) Total number of Samples (stg.sansamples)
\(r =\) Index of the Samples
\(i =\) Index of each parameter
rst.SA.STI - Total order effects
\(S_{Ti}=\frac{V(Y)-V_{Θ_{i}}(E_{Θ_{i}}(Y|Θ_{i}))}{V(Y)}=1-\frac{U_{-i}-E^2(Y)}{V_T(Y)}\)
\(U_{-i}=\frac{1}{n-1}\sum_{r=1}^nf(M_2^r)f(N_i^r)\)
\(E^2(Y)=\frac{1}{n}\sum_{r=1}^nf(M_1^r)f(M_2^r)\)
\(V_T(Y) = \frac{1}{n-1}f^2(M_2^r)-E^2(Y)\)
\(V =\) Variance
\(E(... |...) =\) Conditional expected value
\(Θ =\) Parameters of the model
\(Y =\) Scalar output from the model
\(n =\) Total number of Samples (stg.sansamples)
\(r =\) Index of the Samples
\(i =\) Index of each parameter
References
Model files and folders
“Model Folder name”/ (u-gen)
Folder containing all the model files the model, we recomend using the same name as the repository name for the models we provide but there is no restrictions.
“Source_sbtab_name”.xlsx (u-gen)
Contains the SBtab in .xlsx format.
Matlab/ (u-gen)
Folder containing all model files related to MATLAB®, the model is used in other softwares so here reside all the files that MATLAB® uses or generates
Data/ (a-gen)
Folder containing the model in severall diferent formats relevant for the analysis and files contaning model metadata such as experimental inputs and outputs
model_”model name”.sbproj (a-gen)
Contains the model derived from the SBtab in .sbproj (MATLAB® SimBiology®) format.
model_”model name”.mat (a-gen)
Contains the model derived from the SBtab in .mat format.
data_”model name”.mat (a-gen)
Contains data derived from the SBtab in a .mat format. This data is used to run the model taking into account all the inputs and outputs of the model.
model_”model name”.xml (a-gen)
Contains the model derived from the SBtab in .xml (SBML) format.
Input_”model name”.mat (a-gen)
Contains input data derived from the SBtab in a .mat format for all the experimental inputs.
SBtab_”model name”.mat (a-gen)
Contains the SBtab in .mat format.
Exp/ (a-gen)
Contains a version of the model for each experiment contained in the SBtab. They include all the neccessary inputs and outputs to simulate the supplied experimental conditions.
Model_”model name”_i.mat (a-gen)
Tailor made for the main run of the simulation.
Model_eq_”model name”_i.mat (a-gen)
Tailor made for the equilibration step of the simulation.
Model_detail_”model name”_i.mat (a-gen)
Tailor made for the main run of the simulation. The step size is reduced to generate better graphs
Input_functions/ (a-gen)
Folder containing the functions that are used at run time to give the correct input to all experiments
“model name”_inputi_Ligand.mat (a-gen)
These functions interpolate the input that is supposed to be given to the model at run time.
“model name”_input_creator.mat (a-gen)
Creates the functions above for all experimental inputs.
Results/ (a-gen)
Folder containing the results of the various possible Matlab analysis provided by our workflow
“Analysis name”/ (a-gen)
Each analysis output is stored in its own folder, depending on the analysis run the the results can be saved in either a “Diagnostics”, “Optimization” or “Sensitivity Analysis”. An “Examples” folder is also provided with analysis that were pre-run by us.
“date”/ (a-gen)
The date and time of when the analysis was run, this is auto generated when an user choses to run any alysis.
All_figures.fig (a-gen)
All the plots generated by the analysis stored in a Matlab figure assembly
Analysis.mat (a-gen)
All the data used as input to the analysis, saved as the SBtab and the setting fileconverted to a matlab structs called “sb” and “stg” respectively. And all the outputs generated by running the analysis saved also in a matlab struct called “rst”
“Figure name”.png (a-gen)
All the plots generated by the analysis stored individually as images
Settings/ (u-gen)
Folder containing the settings file
“Settings file name” (u-gen)
Settings file of the model. A place for the user to define all the relevant properties of model simulation that are not stored in SBtab. These are usually things that need to change during optimizations or model development.
tsv/ (a-gen)
Folder containing the SBtab converted to .tsv files for each SBtab that as been run through one of our analysis.
“model name” (a-gen)
Contains the SBtab in .tsv format.
Description of the general terms:
NEURON
Here we describe the usage of the MOD files for two of the provided examples (the third one does not use MOD files).
Nair et al. 2016 (D1 MSN subcellular cascade)
The model requires input in the form of dopamine and calcium. These need to be specified by the user:
Dopamine is set to be 20 nM in assign_calculated_values(). This line should be replaced to whatever the user needs for input. For example, it could be replaced with an expression for dopamine such as the one provided in this mod file, which makes a dopamine pulse with a double exponential shape. The line
DA = 20
should instead read
DA = DA_expression.
The same goes for the calcium input. Alternatively, calcium could be provided via the intracellular concentration of a calcium ion. For example, in this MOD file we use the intracellular calcium concentration due to influx from NMDA receptors, ca_nmdai, which is calculated through a mechanism for calcium accumulation provided in a separate MOD file. Access to the ionic concentrations is provided by NEURON’s “USEION” statement. In this case the variable ca_nmdai needs to be added to the ASSIGNED block.
———IMPORTANT NOTE———
When specifying the calcium input like this care needs to be taken to make sure the units used by NEURON and the units in the model exported to the MOD file match. As mentioned in the article, the models exported as MOD files could have different units for the parameter values from the default units used by NEURON. For instance, NEURON’s default units for internal calculations of ionic concentrations are in millimolars (mM), but the model’s parameters are expressed in nanomolars (nM). It is absolutely paramount to match units, i.e. use the correct scaling for, in this case, the variable ca_nmdai, to provide the model with the right quantity of calcium so that it runs properly:
Ca = ca_nmdai * (1e6)
Viswan et al. 2016 (EGF-stimulated MAPK cascade)
In this model we only use EGF as an input (only this input is used in Figure 7 in Viswan et al. 2018 which we reproduce; otherwise the model may have various other inputs, see the original paper for details [2]).
The input is a step in the concentration of EGF, given by an expression for a sharp sigmoidal function for EGF in assign_calculated_values():
EGF = EGF_level/(1+exp(-(t-EGF_start) * EGF_steepness))
The SBtab format of this model expresses the parameter values in seconds, whereas NEURON’s default unit for time is milliseconds. The conversion script SBtab_to_vfgen provides automatic scaling of time units in SBtab to milliseconds. The concentration units are given in micromolar, and if some species needs to be coupled to a NEURON variable which is expressed in other units (such as NEURON’s default millimolar units), the species or the NEURON variable need to be rescaled as in the example above.
References
Subcellular application
Subcellular application (https://subcellular.humanbrainproject.eu/model/meta) provides a web interface for simulation of biomolecular networks expressed on bionetgen language (https://bionetgen.org/) using network free solver NFsim and reaction-diffusion stochastic systems solver STEPS (http://steps.sourceforge.net/STEPS/documentation.php) Models can be imported from an sbml file. In this repository we used two model examples to exemplify the usage of this tool.
BioNetGen translation of SBtab Nair_2016 model
The model was translated from SBtab model format to rule-based BioNetGen language for the simulation with STEPS and NFsim solvers embedded in the subcellular web app and with the RuleBender
Conversion steps
Run convert_Nair_2016_from_SBTAB_to_SBML.R in RStudio to translate SBtab model to SBML. This step depends on SBtab to SBML converter
Run convert_Nair_2016_from_SBML_to_BNGL.ipynb jupyter notebook to translate from SBML to BioNetGen language
Import the resulted BioNetGen model Nair_2016_optimized_alternative.bngl to the subcellular web app. Add spine geometry .json, .node, .ele, .face files and stimulation pattern stim_DA_complex.tsv. See the subcellular web app help for details
Simulate the final model Nair_2016_optimized_alternative.ebngl in the subcellular web app using STEPS or NFsim solvers
Files and folders
Nair_2016_optimized_alternative.ebngl - extended BNGL model corresponding to the optimized Nair 2016 SBtab model with added geometry and stimulation patterns. Can be imported and simulated in the subcellular web app
SBTAB_Nair_2016 - the folder with the optimized Nair 2016 SBtab model tsv tables.
Nair_2016_optimized.xml - SBML model translated from the optimized Nair 2016 model by convert_Nair_2016_from_SBTAB_to_SBML.R script based on SBtab to SBML converter
Nair_2016_optimized_alternative.bngl - BioNetGen model obtained from Nair_2016_optimized.xml by convert_Nair_2016_from_SBML_to_BNGL.ipynb jupyter notebook which is based on sbml_to_bngl.py conversion tool
spine.ele, spine.face, spine.node, spine.json - these files specify TetGen meshes and model geometry needed for the subcellular web app Geometry section and STEPS solver (see the subcellular web app help for details)
stim_DA_complex.tsv, stim_noDA_complex.tsv - these files specify the stimulation pattern in Simulations section of the subcellular web app (corresponds to the experiments E0 - E9 of the SBtab model).
SBtabVFGEN-master - the folder containing a copy of SBtab to SBML converter
sbml_to_bngl.py - the python tool for conversion of SBML models to BioNetGen language.
BioNetGen translation of SBtab Viswan_2018 model
The model was translated from SBtab model format to rule-based BioNetGen language for the simulation with STEPS and NFsim solvers embedded in the subcellular web app and with the RuleBender
Conversion steps
Run convert_Viswan_2018_for_STEPS_optimised_from_SBTAB_to_SBML.R in RStudio to translate SBtab model to SBML. This step depends on SBtab to SBML converter and LibSBML
Run convert_Viswan_2018_for_STEPS_optimised_from_SBML_to_BNGL.ipynb jupyter notebook to translate from SBML to BioNetGen language. This step depends on sbml_to_bngl.py and on LibSBML or pySB
Import the BioNetGen model (SBTAB_Viswan_2018_alternative.bngl) to the subcellular web app. Add spine geometry ( .json, .node, .ele, .face files) and stimulation pattern (stim_E0.tsv). See the subcellular web app help for details
Simulate final model (SBTAB_Viswan_2018_alternative.ebngl) in the subcellular web app using STEPS or NFsim solvers
Simulate the BioNetGen model with the RuleBender
Files and folders
SBTAB_Viswan_2018_alternative.ebngl - extended BNGL model corresponding to the Viswan_2018_for_STEPS_optimised.xlsx SBTAB model with added geometry and stimulation patterns. Can be imported and simulated in the subcellular web app
Viswan_2018_for_STEPS.xlsx - SBtab model equivalent to the original Viswan_2018 model. It was obtained from Viswan_2018.xlsx by the modification of the model features incompatible with BNGL.
Viswan_2018_for_STEPS_optimised.xlsx - SBtab model equivalent to the optimized Viswan_2018 model. It was obtained from Viswan_2018_optimized.xlsx by the modification of the model features incompatible with BNGL.
SBTAB_Viswan_2018.xml - SBML model translated from Viswan_2018_for_STEPS_optimised.xlsx model by convert_Viswan_2018_for_STEPS_optimised_from_SBTAB_to_SBML.R script based on SBtab to SBML converter
SBTAB_Viswan_2018_alternative.bngl - BioNetGen model obtained from SBTAB_Viswan_2018.xml by convert_Viswan_2018_for_STEPS_optimised_from_SBML_to_BNGL.ipynb jupyter notebook which is based on sbml_to_bngl.py conversion tool
cell.ele, cell.face, cell.node, cell.json - these files specify TetGen meshes and model geometry needed for the subcellular web app Geometry section and STEPS solver (see the subcellular web app help for details)
stim_E0.tsv, stim_E1.tsv - these files specify the stimulation pattern in Simulations section of the subcellular web app (corresponds to the experiments E0 and E1 of the SBtab model).
Viswan_2018_alternative_RuleBender.bngl - BNGL model corresponding to the Viswan_2018_for_STEPS_optimised.xlsx SBtab model with additional section specifying stimulation and BioNetGen solver. Can be imported and simulated in the RuleBender
SBtabVFGEN-master - the folder containing copy of SBtab to SBML converter
SBTAB_Viswan_2018_for_STEPS_optimised - the folder containing tsv tables of Viswan_2018_for_STEPS_optimised.xlsx
sbml_to_bngl.py - the python tool for conversion of SBML models to BioNetGen language.
Conversion of SBML to BioNetGen language
The conversion is implemented in sbml_to_bngl.py python module. Two approaches are supported by sbml_to_bngl.transform() function:
if converter=’pysb’ - the converter based on the Atomizer implemented in pysb.importers.sbml.sbml_translator() function within pySB framework will be used. The Atomizer will try to modify the set of model molecules and reactions to convert them from reaction network to rule-based BioNetGen format.
if converter=’plain’ - a libsbml based converter for sbml level 2, version 4 will be used. This converter produces a bngl approximation to reaction network format of a model. It is assumed that sbml models were obtained by exporting a MATLAB simbiology model to sbml, or by translation of SBTAB model by SBtab to SBML converter.
Models expressed by SBML and SBTAB often are not fully compatible with BNGL. Additional model adaptation steps are required in this case to obtain a working BNGL model. These steps will be partially automatized by sbml_to_bngl.transform() function if adapt_steps argument dictionary ‘list_of_steps’ is nonempty.
The adaptation steps include:
The STEPS and NFsim solvers require different units for species quantities an kinetic rates. An adapted BNGL model provides modifed expressions for all species concentrations and kinetic rates and provides an easy way for units changing by specification of auxilary bngl model parameters: Na and V_comparment_name. These parameters should be selected to: Na=6.022e23 and V_comparment_name = volume of corresponding compartment in liters for NFsim and to: Na=1 and V_comparment_name=1 for STEPS.
Species with fixed concentrations are not supported by NFsim solver. The BNGL model adaptation will modify model reactions such that a fixed species concentration became a model parameter. This parameter can be used for the clamping of species concentration or for the stimulation pattern application
If SBML to BNGL converter implemented in Atomizer is selected then additional transformation steps include renaming of duplicated molecule sites and reparing incorrect molecule names and kinetic rate transformations
In case when MATLAB simbiology is used for SBML model creation, adapt_steps will repare incorrect molecule and parameter names
Compartmental model of BNGL assumes tree structure of the set of model comartments. This assumption is often incompatible with mesh geometries supported by STEPS. The model adaption steps can produce bngl models with flexible compartmental structure
There is a number of incompatibilities between SBTAB and BNGL which still require manual correction. These include concentrations fixed to an expression, functional expressions for reaction rates in case of STEPS solver, some nonstandard types of reaction kinetic functions etc. The adapt_steps detects the cases of known incompatibilities and produces corresponding warning messages
Conversion tools
For this workflow we have developed some conversion tools to facilitate model developement.
The MATLAB® code takes the SBtab file in excel format and generates tab separeted file (.tsv) of this SBtab, an SBML file (.xml) of the model, and two MATLAB® versions of the model (.m and .sbproj). This conversion happens as a setup step of running any of the Matlab analysis, it might be added as a standalone option in the future.
SBtab (.xls,.xlsx) -> Matlab® model (.m .sbproj), SBtab (.tsv), Matlab® SBML (.xml)
We also have R code to perform other conversions in two external repositories:
Code for fixing the SBML produced by MATLAB®
https://github.com/a-kramer/simbiology-sbml-fix
Matlab® SBML (.xml) -> SBML (.xml)
A standalone SBtab to VFgen SBML and NEURON tool
https://github.com/a-kramer/SBtabVFGEN
SBtab (.tsv or .ods) -> VFGEN (.vf) + SBML (.xml) + Neuron (.mod)
Models
We have used multiple models to validate our software tools. Each model has its own repository. In these repositories you can find;
The model in different formats relevant for the various tools that we provide
“Matlab” folder with:
Settings file relevant to use the model with our MATLAB® tools
Some examples of the output provided by our MATLAB® tools after running an Analysis
This is also the folder were all the run time ouptuts of the matlab analysis are stored.
“Bionetgen and Steps” folder with relevant files for running the model in the subcellular aplication
“Neuron” folder with relevant files for running the model in NEURON.
Note: Model_Fujita_2010 does not contain “Bionetgen and Steps” and “Neuron” folders as this model was not implemented in these softwares.
Links to the model repositories:
Fujita et al. 2010
Model by Fujita et al 2010 as one of the proposed benchmark models provided by Hass et al. 2019. The model represents epidermal growth factor (EGF)-dependent Akt pathway. Data from [Fig. 1b] with a corresponding EGF step input from [Fig. 1a] in Fujita et al. 2010 was used for estimation and Global Sensitity Analysis of 12 parameters. We performed parameter estimation starting from the model parameters provided in the publication, using a search space 5 order of magnitude above and below for each parameter.
Reactions chosen for parameter estimation and Global Sensitivity Analysis
EGFR + EGF <=> EGF_EGFR
pEGFR + Akt <=> pEGFR_Akt
pEGFR_Akt -> pEGFR + pAkt
pEGFR -> null
pAkt + S6 <=> pAkt_S6
pAkt_S6 -> pAkt + pS6
pAkt -> Akt
pS6 -> S6
EGF_EGFR -> pEGFR
EGFR -> null
Model Repository
https://github.com/jpgsantos/Model_Fujita_2010
References
Nair et al. 2016
We illustrate the Subcellular Workflow with a model depicting the emergence of the eligibility trace observed in reinforcement learning in striatal D1 medium spiny neurons (D1 MSN) (Nair et al. 2016). Here, an intracellular increase in calcium representing excitatory synaptic input leads to synaptic potentiation only when it is followed by a reinforcing dopamine input. These two signaling cascades, starting with a calcium train and a dopamine transient, are illustrated in Fig. 1a. The first pathway (depicted in blue) represents calcium-dependent activation of Ca2+/calmodulin-dependent protein kinase II (CaMKII), its autophophorylation, and the phosphorylation of a generic CaMKII substrate that represents long term potentiation (LTP). In the second pathway (species in red), dopamine initiates a G-protein dependent cascade which results in the phosphorylation of dopamine- and cAMP-regulated phosphoprotein, 32 kDa (DARPP-32) turning into an inhibitor of protein phosphatase 1 (PP1) that otherwise dephosphorylates CaMKII and its substrate. Substrate phosphorylation is maximal when the time window between the calcium and dopamine inputs is sufficiently short (input-interval constraint mediated by DARPP-32 via PP1 inhibition), and when intracellular calcium elevation is followed by dopamine (input-order constraint mediated by another phosphoprotein, the cyclic AMP-regulated phosphoprotein, 21 kDa (ARPP-21) that sequesters calcium/calmodulin if dopamine arrives before calcium).
CaMKII is autophosphorylated both in the cytosol and the post synaptic density (PSD) with a custom MATLAB rate function as described in Li et al. (2012). To run the model in different software, we substitute the autophosphorylation rate function with the same set of bimolecular reactions (simplified version of reactions from Pepke et al. 2010) for both compartments. The reactions represent the formation of a complex with two fully activated CaMKII monomers and a catalytic step in which one monomer phosphorylates the other. Schematics can be found in Fig. 1b along with the six new parameters. We estimated the parameters using simulated data from the published model with simulation setups in Fig. 1c depicting different timings of the dopamine input relative to the calcium input.
Figure 1. (A) Simplified schematics of the model. Species in the calcium cascade are depicted in blue, and species in the dopamine cascade are depicted red. Simulated time course data of the species with a beige background were used in parameter estimation. Figure adapted from Nair et al. (2016). (B) Illustration of the bimolecular reactions in CaMKII autophosphorylation with yellow circles representing activated CaMKII, and blue circles representing phosphate groups. Newly introduced parameters with their IDs in the updated model are shown above/below the arrows. (C) Timing of the dopamine input (Δt ={-4,-3,-2,-1,0,1,2,3,4} corresponding to experiments E0-E8) relative to calcium increase (time 0), and a single experiment without dopamine (E9).
Reactions chosen for parameter estimation and Global Sensitivity Analysis
Four reactions representing autophosphorylation in one compartment.
Reactions in cytosol:
CaMKII_CaM_Ca4 + CaMKII_CaM_Ca4 <=> CaMKII_CaM_Ca4_CaMKII_CaM_Ca4
CaMKII_CaM_Ca4_CaMKII_CaM_Ca4 -> pCaMKII_CaM_Ca4 + CaMKII_CaM_Ca4
pCaMKII_CaM_Ca4 + CaMKII_CaM_Ca4 <=> pCaMKII_CaM_Ca4_CaMKII_CaM_Ca4
pCaMKII_CaM_Ca4_CaMKII_CaM_Ca4 -> pCaMKII_CaM_Ca4 + pCaMKII_CaM_Ca4
Reactions in PSD:
CaMKII_CaM_Ca4_psd + CaMKII_CaM_Ca4_psd <=> CaMKII_CaM_Ca4_psd_CaMKII_CaM_Ca4_psd
CaMKII_CaM_Ca4_psd_CaMKII_CaM_Ca4_psd -> pCaMKII_CaM_Ca4_psd + CaMKII_CaM_Ca4_psd
pCaMKII_CaM_Ca4_psd + CaMKII_CaM_Ca4_psd <=> pCaMKII_CaM_Ca4_psd_CaMKII_CaM_Ca4_psd
pCaMKII_CaM_Ca4_psd_CaMKII_CaM_Ca4_psd -> pCaMKII_CaM_Ca4_psd + pCaMKII_CaM_Ca4_psd
Model Repository
https://github.com/jpgsantos/Model_Nair_2016
References
Viswan et al. 2018
Model from the FindSim framework by Viswan et al. 2018. The model represents an epidermal growth factor (EGF)-dependent mitogen-activated protein kinase (MAPK) signaling pathway (the green block from [Fig. 7a]) which measures MAPK phosphorylation. Here, two simulation experiments with EGF step inputs of different sizes are used to perform parameter estimation on 29 model parameters corresponding to reactions involved in MAPK phosphorylation (see below). For this we used activated MAPK curves in [Fig. 7b and 7c].
Reactions chosen for parameter estimation and Global Sensitivity Analysis
GTP_Ras + craf_1_p <=> Raf_p_GTP_Ras
GEF_p -> inact_GEF
GTP_Ras -> GDP_Ras
GAP_p -> GAP
MAPK_p_p + craf_1_p <=> MAPK_p_p_feedback_cplx
MAPK_p_p_feedback_cplx -> MAPK_p_p + craf_1_p_p
MAPKK_p_p + MAPK <=> MAPKKtyr_cplx
MAPKKtyr_cplx -> MAPKK_p_p + MAPK_p
MAPKK_p_p + MAPK_p <=> MAPKKthr_cplx
MAPKKthr_cplx -> MAPKK_p_p + MAPK_p_p
MAPKK + Raf_p_GTP_Ras <=> Raf_p_GTP_Ras_1_cplx
Raf_p_GTP_Ras_1_cplx -> MAPKK_p + Raf_p_GTP_Ras
MAPKK_p + Raf_p_GTP_Ras <=> Raf_p_GTP_Ras_2_cplx
Raf_p_GTP_Ras_2_cplx -> MAPKK_p_p + Raf_p_GTP_Ras
inact_GEF + GDP_Ras <=> basal_GEF_activity_cplx
basal_GEF_activity_cplx -> inact_GEF + GTP_Ras
GEF_p + GDP_Ras <=> GEF_p_act_Ras_cplx
GEF_p_act_Ras_cplx -> GEF_p + GTP_Ras
GAP + GTP_Ras <=> GAP_inact_Ras_cplx
GAP_inact_Ras_cplx -> GAP + GDP_Ras
Model Repository
https://github.com/jpgsantos/Model_Viswan_2018