Simulation and Scoring

f_sim_score

Code

function [score,rst,rst_not_simd] = f_sim_score(parameters,stg,mmf)

%Turn off Dimension analysis warning from simbiology
warning('off','SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless')

% Call the function that simulates the model
rst = f_prep_sim(parameters,stg,mmf);

% Call the function that scores
rst = f_score(rst,stg,mmf);

% Get the total score explicitly for optimization functions
score = rst.st;

rst_not_simd = rmfield( rst , 'simd');
end

Calls the function that runs the simulations and the function that scores the output of the runs.

• Inputs

  • stg

  • parameters - (double) Set of parameters that we are working on

• Outputs

  • score - rst.st

  • rst - rst.simd, rst.xfinal, rst.sd, rst.se, and rst.st

  • rst_not_simd - rst.xfinal, rst.sd, rst.se, and rst.st

• Calls - f_prep_sim, f_score

• Loads

f_prep_sim

Code

function rst = f_prep_sim(parameters,stg,mmf)

% Save variables that need to be mantained over multiple function calls
% persistent modelobj
persistent sbtab
persistent Data

data_model = mmf.model.data.data_model;

% Import the data on the first run
if isempty(sbtab)
    
    %Find correct path for loading depending on the platform
    load(data_model,'Data','sbtab')
end

% Set the parameters that are going to be used for the simulation to the
% default ones as definded in the SBTAB
rt.par(:,1) = [sbtab.defpar{:,2}];

% Check if the parametrer needs to be set to the value relevant for Profile
% Likelihood
if isfield(stg,"PLind") && isfield(stg,"PLval")
    parameters = [parameters(1:stg.PLind-1) stg.PLval parameters(stg.PLind:end)];
end

% Iterate over all the parameters of the model
for n = 1:size(rt.par,1)
    
    % Check that a parameter should be changed from default
    if stg.partest(n) > 0
        
        % Set the parameters are being tested
        rt.par(n) = 10.^(parameters(stg.partest(n,1)));
    end
    
    if isfield(stg,'tci')
        
        % Check that there are thermodynamic constraints to implement
        if ~isempty(stg.tci)
            
            % Choose the parameters that need to be calculated with other
            % parameters due to thermodynamic constraints
            if ismember(n,stg.tci)
                
                % Check that a parameter should be changed from default
                if stg.partest(n) > 0
                    
                    % Iterate over the parameters that need to be mutiplied
                    % for calculating the parameter that depends on the
                    % thermodynamic constraints
                    for m = 1:size(stg.tcm,2)
                        
                        % Check that the parameter that is going to be used
                        % to calculate the parameter dependent on
                        % thermodynamic constraintsis is not the default
                        if stg.partest(stg.tcm(n,m),1) > 0
                            
                            % Check if the parametrer needs to be set to
                            % the value relevant for Profile Likelihood
                            if isfield(stg,"PLind")
                                if stg.partest(stg.tcm(n,m),1) ==...
                                        stg.PLind
                                    parameters(stg.partest(...
                                        stg.tcm(n,m),1))...
                                        = stg.PLval;
                                end
                            end
                            
                            % Make the appropriate multiplications to get
                            % the thermodinamicly constrained parameter
                            rt.par(n) = rt.par(n).*(10.^...
                                (parameters(stg.partest(...
                                stg.tcm(n,m),1))));
                        else
                            
                            % Make the appropriate multiplications to get
                            % the thermodinamicly constrained parameter
                            rt.par(n) = rt.par(n).*...
                                (sbtab.defpar{stg.tcm(n,m),2});
                        end
                    end
                    
                    % Iterate over the parameters that need to be divided
                    % for calculating the parameter that depends on the
                    % thermodynamic constraints
                    for m = 1:size(stg.tcd,2)
                        
                        % Check that the parameter that is going to be used
                        % to calculate the parameter dependent on
                        % thermodynamic constraintsis is not the default
                        if stg.partest(stg.tcd(n,m),1) > 0
                            
                            % Check if the parametrer needs to be set to
                            % the value relevant for Profile Likelihood
                            if isfield(stg,"PLind")
                                if stg.partest(stg.tcd(n,m),1) ==...
                                        stg.PLind
                                    parameters(stg.partest(...
                                        stg.tcd(n,m),1))...
                                        = stg.PLval;
                                end
                            end
                            % Make the appropriate divisions to get the
                            % thermodinamicly constrained parameter
                            rt.par(n) = rt.par(n)./(10.^...
                                (parameters(stg.partest(...
                                stg.tcd(n,m),1))));
                        else
                            
                            % Make the appropriate divisions to get the
                            % thermodinamicly constrained parameter
                            rt.par(n) = rt.par(n)./...
                                (sbtab.defpar{stg.tcd(n,m),2});
                        end
                    end
                end
            end
        end
    end
end

%  Set the start amount for the species in the model to 0
rt.ssa = zeros(size(sbtab.species,1),max(stg.exprun));

% Initialize the results variable
rst = [];
rst.parameters = rt.par;
% Iterate over all the experiments that are being run
for n = stg.exprun
    
    % Try catch used because iterations errors can happen unexectedly and
    % we want to be able to continue simulations
    try
        
        % If the correct setting is chosen display messages to console
        if stg.simcsl
            disp("Running dataset number " + n +...
                " of " + stg.exprun(end))
        end
        
        % Check that this is not the first time the experiments are being
        % run and that the start values for the species are different from
        % the previous experiment
        if n ~= stg.exprun(1) && ...
                min([sbtab.datasets(n).start_amount{:,2}] ==...
                [sbtab.datasets(max(1,stg.exprun(...
                find(stg.exprun==n)-1))).start_amount{:,2}])
            
            
            % Set the values of the start amounts to the values obtained
            % after the first equilibration
            rt.ssa(:,n) =...
                rt.ssa(:,stg.exprun(find(stg.exprun==n)-1));
            if stg.simdetail
                rt.ssa(:,n+2*stg.expn) =...
                    rt.ssa(:,stg.exprun(find(stg.exprun==n)-1));
            end
        else
            
            % Iterate over the numbre of species that need a starting value
            % different than 0
            for j = 1:size(sbtab.datasets(n).start_amount,1)
                
                % Set the start amount of the species to the number defined
                % in the sbtab for each experiment
                rt.ssa(sbtab.datasets(n...
                    ).start_amount{j,3},n+stg.expn) =...
                    sbtab.datasets(n).start_amount{j,2};
            end
            
            % Equilibrate the model
            rst = f_sim(n+stg.expn,stg,rt,rst,mmf);
            
            for j = 1:size(sbtab.species,1)
                
                % Set the starting amount for species that after
                % equilibrium have very low values to zero
                if rst.simd{n+stg.expn}.Data(end,j) < 1.0e-15
                    rt.ssa(j,n) = 0;
                    
                    % Set the starting amount for the rest of the species
                else
                    rt.ssa(j,n) =...
                        rst.simd{n+stg.expn}.Data(end,j);
                    if stg.simdetail
                        rt.ssa(j,n+2*stg.expn) =...
                            rst.simd{n+stg.expn}.Data(end,j);
                    end
                end
            end
        end
        
        % Simulate the model
        rst = f_sim(n,stg,rt,rst,mmf);
        
        try
            if stg.simdetail
                rst = f_sim(n+2*stg.expn,stg,rt,rst,mmf);
            end
        catch
        end
        
        % Check If the times of the simultaion output and the simulation
        % data from SBTAB match, if they don't it means that the simulator
        % didn't had enough time to run the model (happens in some
        % unfavorable configurations of parameters, controlled by stg.maxt
        if size(Data(n).Experiment.t,1) ~=...
                size(rst.simd{n}.Data(:,end),1)
            
            % Set the simulation output to be 0, this is a non function
            % value that the score function expects in simulations that did
            % not worked properly
            rst.simd{n} = 0;
        end
    catch
        
        % Set the simulation output to be 0, this is a non function value
        % that the score function expects in simulations that did not
        % worked properly
        rst.simd{n} = 0;
    end
end
end

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

  • rst - rst.simd

• Calls - f_sim

• Loads - data.mat, model.mat

f_sim

Code

function rst = f_sim(exp_n,stg,rt,rst,mmf)

% Save variables that need to be mantained over multiple function calls
persistent model_run
persistent config_run

% If the function is called for the first time, load the appropriate model
% and compile the code for simulation run
if isempty(model_run)
    
    warning('off','SimBiology:InvalidSpeciesInitAmtUnits')
    
    %Generate an empty array to be populated with the model suited for each
    %equilibration and experiment%
    model_run = cell(1,stg.expn*2);
    config_run = cell(1,stg.expn*2);
    
    model_exp_default = mmf.model.data.model_exp.default;
    model_exp_eq = mmf.model.data.model_exp.equilibration;
    model_exp_detail = mmf.model.data.model_exp.detail;
    
    % Iterate over the experiments thar are being run
    for n = stg.exprun
        
        if stg.simdetail
            % Load the models for equilibrium
            load(model_exp_detail + n + ".mat",'model_exp','config_exp')

            % Place the loaded models in the correct place in the array,
            % models for equilibrium are set to be on the second half of
            % the array
            model_run{n+2*stg.expn} = model_exp;
            config_run{n+2*stg.expn} = config_exp;
            
            % Compile the matlab code that is going to simulate the model
            % using matlab built in function if the option is selected in
            % settings
            if stg.sbioacc
                sbioaccelerate(model_run{n+2*stg.expn},config_run{n+2*stg.expn});
            end
        end
        % Load the models for equilibrium
        load(model_exp_eq + n + ".mat",'model_exp','config_exp')

        % Place the loaded models in the correct place in the array, models
        % for equilibrium are set to be on the second half of the array
        model_run{n+stg.expn} = model_exp;
        config_run{n+stg.expn} = config_exp;
        
        % Compile the matlab code that is going to simulate the model using
        % matlab built in function if the option is selected in settings
        if stg.sbioacc
            sbioaccelerate(model_run{n+stg.expn},config_run{n+stg.expn});
        end
        
        % Load the models for main run
        load(model_exp_default + n + ".mat",'model_exp','config_exp')
        
        % Place the loaded models in the correct place in the array, models
        % for main run are set to be on the first half of the array
        model_run{n} = model_exp;
        config_run{n} = config_exp;
        
        % Compile the matlab code that is going to simulate the model using
        % matlab built in function if the option is selected in settings
        if stg.sbioacc
            sbioaccelerate(model_run{n},config_run{n});
        end
    end
end

% substitute the start amount of the species in the model with the correct
% ones for  simulations
set(model_run{exp_n}.species(1:size(rt.ssa(:,exp_n),1)),{'InitialAmount'},num2cell(rt.ssa(:,exp_n)))

% Substitute the values of the parameters in the model for the correct one
% for simultaions
set(model_run{exp_n}.parameters(1:size(rt.par,1)),{'Value'},num2cell(rt.par))

%simulate the model using matlab built in function
rst.simd{exp_n} = sbiosimulate(model_run{exp_n},config_run{exp_n});
end

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.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.

  • rst - rst.simd

• Outputs

  • rst - rst.simd

• Calls - Sbioaccelerate, Sbiosimulate

• Loads - Ready to run model, Ready to run model equilibration

f_score

Code

function rst = f_score(rst,stg,mmf)

persistent sbtab
persistent Data

data_model = mmf.model.data.data_model;

% Import the data on the first run
if isempty(Data)
    load(data_model,'Data','sbtab')
end

% Iterate over the number of experiments
for n = stg.exprun
    
    % Iterate over the number of datasets per experiment
    for j = 1:size(sbtab.datasets(n).output,2)
        
        % Calculate score per dataset if there are no errors
        if rst.simd{n} ~= 0
            
            data = Data(n).Experiment.x(:,j);
            data_sd = Data(n).Experiment.x_SD(:,j);
            number_points = size(Data(n).Experiment.x(:,j),1);
            sim_results = f_normalize(rst,stg,n,j,mmf);
            rst.xfinal{n,1}(j) = sim_results(end);
            
            % Calculate score using formula that accounts for normalization
            % with the starting point of the result
            if stg.useLog == 4
                rst.sd(j,n) = sum(((data-sim_results)./...
                    (data_sd*sqrt(number_points))).^2);
                 rst.sdtest(j,n) = sum(((data-sim_results)./...
                     (data_sd*sqrt(number_points))).^2);

%             elseif stg.useLog == 5  
%                 rst.sd{n,1}(j) = sum(((data-sim_results)./...
%                         (data_sd)).^2);
            else
                    rst.sd(j,n) = sum(((data-sim_results)./...
                        (data_sd)).^2)/(number_points);
                    rst.sdtest(j,n) = sum(((data-sim_results)./...
                        (data_sd)).^2)/(number_points);
            end
            
            % If there are errors output a very high score value (10^10)
        elseif rst.simd{n} == 0 || rst.sd(n,j) == inf
            
                rst.sd(n,j) = stg.errorscore;
                rst.xfinal{n,1}(j) = 0;
%                 rst.sdtest(j,n) = stg.errorscore;
        end
%         rst.sd{n,1}(j)
%         rst.sdtest(j,n)
% max(0,log10(rst.sd{n,1}(j)));
% max(0,log10(rst.sdtest(j,n)));

        % Calculate the log10 of dataset score if option selected
        if stg.useLog == 1
            rst.sd(j,n) = max(0,log10(rst.sd(j,n)));
        end
    end

    % Calculate score per experiment
    rst.se(n,1) = sum(rst.sd(:,n));
%     rst.setest(n,1) = sum(rst.sdtest(:,n));
    
    % Calculate the log10 of experiment score if option selected
    if stg.useLog == 2
        rst.se(n,1) = log10(rst.se(n,1));
    end
end

% Calculate score per experiment
rst.st = sum(rst.se);

% Calculate the log10 of total score if option selected
if stg.useLog == 3
    rst.st = log10(rst.st);
end
end

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

  • rst - rst.simd

  • stg - stg.name, stg.exprun, stg.useLog

• Outputs

  • rst.st - rst.xfinal, rst.sd, rst.se, rst.st

• Calls

• Loads - data.mat