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.folder_model, 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 %Generate an empty array to be populated with the model suited for each 12 %equilibration and experiment% 13 model_run = cell(1,stg.expn*2); 14 config_run = cell(1,stg.expn*2); 15 16 model_exp_default = mmf.model.data.model_exp.default; 17 model_exp_eq = mmf.model.data.model_exp.equilibration; 18 model_exp_detail = mmf.model.data.model_exp.detail; 19 20 % Iterate over the experiments thar are being run 21 for n = stg.exprun 22 23 if stg.simdetail 24 % Load the models for equilibrium 25 load(model_exp_detail + n + ".mat",'model_exp','config_exp') 26 27 % Place the loaded models in the correct place in the array, 28 % models for equilibrium are set to be on the second half of 29 % the array 30 model_run{n+2*stg.expn} = model_exp; 31 config_run{n+2*stg.expn} = config_exp; 32 33 % Compile the matlab code that is going to simulate the model 34 % using matlab built in function if the option is selected in 35 % settings 36 if stg.sbioacc 37 sbioaccelerate(model_run{n+2*stg.expn},config_run{n+2*stg.expn}); 38 end 39 end 40 % Load the models for equilibrium 41 load(model_exp_eq + n + ".mat",'model_exp','config_exp') 42 43 % Place the loaded models in the correct place in the array, models 44 % for equilibrium are set to be on the second half of the array 45 model_run{n+stg.expn} = model_exp; 46 config_run{n+stg.expn} = config_exp; 47 48 % Compile the matlab code that is going to simulate the model using 49 % matlab built in function if the option is selected in settings 50 if stg.sbioacc 51 sbioaccelerate(model_run{n+stg.expn},config_run{n+stg.expn}); 52 end 53 54 % Load the models for main run 55 load(model_exp_default + n + ".mat",'model_exp','config_exp') 56 57 % Place the loaded models in the correct place in the array, models 58 % for main run are set to be on the first half of the array 59 model_run{n} = model_exp; 60 config_run{n} = config_exp; 61 62 % Compile the matlab code that is going to simulate the model using 63 % matlab built in function if the option is selected in settings 64 if stg.sbioacc 65 sbioaccelerate(model_run{n},config_run{n}); 66 end 67 end 68end 69 70% substitute the start amount of the species in the model with the correct 71% ones for simulations 72set(model_run{exp_n}.species(1:size(rt.ssa(:,exp_n),1)),{'InitialAmount'},num2cell(rt.ssa(:,exp_n))) 73 74% Substitute the values of the parameters in the model for the correct one 75% for simultaions 76set(model_run{exp_n}.parameters(1:size(rt.par,1)),{'Value'},num2cell(rt.par)) 77 78%simulate the model using matlab built in function 79rst.simd{exp_n} = sbiosimulate(model_run{exp_n},config_run{exp_n}); 80end
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)
stg - stg.expn, stg.folder_model, stg.name, stg.sbioacc
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{n,1}(j) = sum(((data-sim_results)./... 32 (data_sd*sqrt(number_points))).^2); 33 else 34 rst.sd{n,1}(j) = sum(((data-sim_results)./... 35 (data_sd)).^2)/(number_points); 36 end 37 38 % If there are errors output a very high score value (10^10) 39 elseif rst.simd{n} == 0 || rst.sd{n,1}(j) == inf 40 41 rst.sd{n,1}(j) = stg.errorscore; 42 rst.xfinal{n,1}(j) = 0; 43 end 44 45 % Calculate the log10 of dataset score if option selected 46 if stg.useLog == 1 47 rst.sd{n,1}(j) = max(0,log10(rst.sd{n,1}(j))); 48 end 49 end 50 51 % Calculate score per experiment 52 rst.se(n,1) = sum(rst.sd{n,1}); 53 54 % Calculate the log10 of experiment score if option selected 55 if stg.useLog == 2 56 rst.se(n,1) = log10(rst.se(n,1)); 57 end 58end 59 60% Calculate score per experiment 61rst.st = sum(rst.se); 62 63% Calculate the log10 of total score if option selected 64if stg.useLog == 3 65 rst.st = log10(rst.st); 66end 67end
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.
Inputs
stg - stg.folder_model, stg.name, stg.exprun, stg.useLog
Outputs
rst.st - rst.xfinal, rst.sd, rst.se, rst.st
Calls
Loads - data.mat