sbiosteadystate

Find steady state of SimBiology model

Syntax

[success, variant_out] = sbiosteadystate(model)

[success, variant_out] = sbiosteadystate(model, variant_in)

[success, variant_out] = sbiosteadystate(model, variant_in, scheduleDose)

[success, variant_out, model_out] = sbiosteadystate(model,___)

[success, variant_out, model_out, exitInfo] = sbiosteadystate(model,___)

[___] = sbiosteadystate(___, Name,Value)

Description

[success, variant_out] = sbiosteadystate(model) attempts to find a steady state of a SimBiology^® model, model. The function returns success, which is true if a steady state was found, and a SimBiology Variant object, variant_out, with all non-constant species, compartments, and parameters of the model having the steady-state values. If a steady state was not found, then the success is false and variant_out contains the last values found by the algorithm.

[success, variant_out] = sbiosteadystate(model, variant_in) applies the alternate quantity values stored in a SimBiology variant object, variant_in, to the model before trying to find the steady-state values.

[success, variant_out] = sbiosteadystate(model, variant_in, scheduleDose) applies a SimBiology schedule dose object, scheduleDose, or a vector of schedule doses to the corresponding model quantities before trying to find the steady state values. Only doses at time = 0 are allowed, that is, the dose time of each dose object must be 0. To specify a dose without specifying a variant, set variant_in to an empty array, [].

[success, variant_out, model_out] = sbiosteadystate(model,___) also returns a SimBiology model, model_out that is a copy of the input model with the states set to the steady-state solution that was found. Also, model_out has all initial assignment rules disabled.

[success, variant_out, model_out, exitInfo] = sbiosteadystate(model,___) also returns the exit information about the steady state computation.

[___] = sbiosteadystate(___, Name,Value) uses additional options specified by one or more Name,Value pair arguments.

Examples

Find a steady state of a simple gene regulation model

This example shows how to find a steady state of a simple gene regulation model, where the protein product from translation controls transcription.

Load the sample SimBiology project containing the model, m1. The model has five reactions and four species.

sbioloadproject('gene_reg.sbproj','m1')

Display the model reactions.

m1.Reactions

ans = 
   SimBiology Reaction Array

   Index:    Reaction:                    
   1         DNA -> DNA + mRNA            
   2         mRNA -> mRNA + protein       
   3         DNA + protein <-> DNA_protein
   4         mRNA -> null                 
   5         protein -> null              

A steady state calculation attempts to find the steady state values of non-constant quantities. To find out which model quantities are non-constant in this model, use sbioselect.

sbioselect(m1,'Where','Constant*','==',false)

ans = 
   SimBiology Species Array

   Index:    Compartment:    Name:          Value:    Units:  
   1         unnamed         DNA            50        molecule
   2         unnamed         DNA_protein    0         molecule
   3         unnamed         mRNA           0         molecule
   4         unnamed         protein        0         molecule

There are four species that are not constant, and the initial amounts of three of them are set to zero.

Use sbiosteadystate to find the steady state values for those non-constant species.

[success,variantOut] = sbiosteadystate(m1)

success = logical
   1

variantOut = 
   SimBiology Variant - SteadyState (inactive)

   ContentIndex:     Type:        Name:             Property:           Value:
   1                 compartment  unnamed           Capacity            1
   2                 species      DNA               InitialAmount       8.79024
   3                 species      DNA_protein       InitialAmount       41.2098
   4                 species      mRNA              InitialAmount       1.17203
   5                 species      protein           InitialAmount       23.4406
   6                 parameter    Transcription.k1  Value               0.2
   7                 parameter    Translation.k2    Value               20
   8                 parameter    [Binding/Unbin... Value               0.2
   9                 parameter    [Binding/Unbin... Value               1
   10                parameter    [mRNA Degradat... Value               1.5
   11                parameter    [Protein Degra... Value               1

The initial amounts of all species of the model have been set to the steady-state values. DNA is a conserved species since the total of DNA and DNA_protein is equal to 50.

You can also use a variant to store alternate initial amounts and use them during the steady state calculation. For instance, you could set the initial amount of DNA to 100 molecules instead of 50.

variantIn = sbiovariant('v1');
addcontent(variantIn,{'species','DNA','InitialAmount',100});
[success2,variantOut2,m2] = sbiosteadystate(m1,variantIn)

success2 = logical
   1

variantOut2 = 
   SimBiology Variant - SteadyState (inactive)

   ContentIndex:     Type:        Name:             Property:           Value:
   1                 compartment  unnamed           Capacity            1
   2                 species      DNA               InitialAmount       12.7876
   3                 species      DNA_protein       InitialAmount       87.2124
   4                 species      mRNA              InitialAmount       1.70502
   5                 species      protein           InitialAmount       34.1003
   6                 parameter    Transcription.k1  Value               0.2
   7                 parameter    Translation.k2    Value               20
   8                 parameter    [Binding/Unbin... Value               0.2
   9                 parameter    [Binding/Unbin... Value               1
   10                parameter    [mRNA Degradat... Value               1.5
   11                parameter    [Protein Degra... Value               1

m2 = 
   SimBiology Model - cell 

   Model Components:
     Compartments:      1
     Events:            0
     Parameters:        6
     Reactions:         5
     Rules:             0
     Species:           4
     Observables:       0

Since the algorithm has found a steady state, the third output m2 is the steady state model, where the values of non-constant quantities have been set to steady state values. In this example, the initial amounts of all four species have been updated to steady state values.

m2.Species

ans = 
   SimBiology Species Array

   Index:    Compartment:    Name:          Value:     Units:  
   1         unnamed         DNA            12.7876    molecule
   2         unnamed         DNA_protein    87.2124    molecule
   3         unnamed         mRNA           1.70502    molecule
   4         unnamed         protein        34.1003    molecule

Input Arguments

model — SimBiology model

SimBiology model object

SimBiology model, specified as a SimBiology Model object.

variant_in — SimBiology variant

variant object

SimBiology variant, specified as a Variant object. The alternate quantity values stored in the variant are applied to the model before finding the steady state.

scheduleDose — Dosing information

SimBiology schedule dose object

Dosing information, specified as a SimBiology schedule dose object. The dose must be bolus, that is, there must be no time lag or administration time for the dose. In other words, its LagParameterName and DurationParameterName properties must be empty, and the dose time (the Time property) must be 0. For details on how to create a bolus dose, see Creating Doses Programmatically.

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'AbsTol',1e-6 specifies to use the absolute tolerance value of 10–6.

'Method' — Method to compute steady state

'auto' (default) | 'simulation' | 'algebraic'

Method to compute the steady state of model, specified as the comma-separated pair consisting of 'Method' and a character vector 'auto', 'simulation', or 'algebraic'. The default ('auto') behavior is to use the 'algebraic' method first. If that method is unsuccessful, the function uses the 'simulation' method.

For the simulation method, the function simulates the model and uses finite differencing to detect a steady state. For details, see Simulation Method.

For the algebraic method, the function computes a steady state by finding a root of the flux function algebraically. For nonlinear models, this method requires Optimization Toolbox™. For details, see Algebraic Method.

Note

The steady state returned by the algebraic method is not guaranteed to be the same as the one found by the simulation method. The algebraic method is faster since it involves no simulation, but the simulation method might be able to find a steady state when the algebraic method could not.

Example: 'Method','algebraic'

'AbsTol' — Absolute tolerance to detect convergence

1e-8 (default) | positive, real scalar

Absolute tolerance to detect convergence, specified as the comma-separated pair consisting of 'AbsTol' and a positive, real scalar.

When you use the algebraic method, the absolute tolerance is used to specify optimization settings and detect convergence. For details, see Algebraic Method.

When you use the simulation method, the absolute tolerance is used to determine convergence when finding a steady state solution by forward integration as follows: $$\left(\Vert \frac{d\overrightarrow{S}}{dt}\Vert <AbsTol\right)or\left(\Vert \frac{d\overrightarrow{S}}{dt}\Vert <RelTol\ast \Vert \overrightarrow{S}\Vert \right)$$, where $$\overrightarrow{S}$$ is a vector of nonconstant species, parameters, and compartments.

'RelTol' — Relative tolerance to detect convergence

1e-6 (default) | positive, real scalar

Relative tolerance to detect convergence, specified as the comma-separated pair consisting of 'RelTol' and a positive, real scalar. This name-value pair argument is used for the simulation method only. The algorithm converges and reports a steady state if the algorithm finds model states by forward integration, such that $$\left(\Vert \frac{d\overrightarrow{S}}{dt}\Vert <AbsTol\right)or\left(\Vert \frac{d\overrightarrow{S}}{dt}\Vert <RelTol\ast \Vert \overrightarrow{S}\Vert \right)$$, where $$\overrightarrow{S}$$ is a vector of non-constant species, parameters, and compartments.

'MaxStopTime' — Maximum amount of simulation time to take before terminating without a steady state

100000 (default) | positive integer

Maximum amount of simulation time to take before terminating without a steady state, specified as the comma-separated pair consisting of 'MaxStopTime' and a positive integer. This name-value pair argument is used for the simulation method only.

'MinStopTime' — Minimum amount of simulation time to take before searching for a steady state

1 (default) | positive integer

Minimum amount of simulation time to take before searching for a steady state, specified as the comma-separated pair consisting of 'MinStopTime' and a positive integer. This name-value pair argument is used for the simulation method only.

Output Arguments

success — Flag to indicate if a steady state of the model is found

true | false

Flag to indicate if a steady state of the model is found, returned as true or false.

variant_out — SimBiology variant

variant object

SimBiology variant, returned as a variant object. The variant includes all species, parameters, and compartments of the model with the non-constant quantities having the steady-state values.

model_out — SimBiology model at the steady state

model object

SimBiology model at the steady state, returned as a model object. model_out is a copy of the input model, with the non-constant species, parameters, and compartments set to the steady-state values. Also, model_out has all initial assignment rules disabled. Simulating the model at steady state requires that initial assignment rules be inactive, since these rules can modify the values in variant_out.

Note

• If you decide to commit the variant_out to the input model that has initial assignment rules, then model is not expected to be at the steady state because the rules perturb the system when you simulate the model.

• model_out is at steady state only if simulated without any doses.

exitInfo — Exit information about steady state computation

character vector

Exit information about the steady state computation, returned as a character vector. The information contains different messages for corresponding exit conditions.

• Steady state found (simulation) – A steady state is found using the simulation method.

• Steady state found (algebraic) – A steady state is found using the algebraic method.

• Steady state found (unstable) – An unstable steady state is found using the algebraic method.

• Steady state found (possibly underdetermined) – A steady state that is, possibly, not asymptotically stable is found using the algebraic method.

• No Steady state found – No steady state is found.

• Optimization Toolbox (TM) is missing – The method is set to 'algebraic' for nonlinear models and Optimization Toolbox is missing.

More About

Simulation Method

sbiosteadystate simulates the model until MaxStopTime. During the simulation, the function approximates the gradient using finite differencing (forward difference) over time to detect a steady state.

Algebraic Method

sbiosteadystate tries to find a steady state of the model algebraically by finding a root of the flux function v. The flux function includes reaction equations, rate rules, and algebraic equations, that is, v(X,P) = 0, where X and P are nonconstant quantities and parameters of the model. Thereby the mass conservation imposed by the reaction equations is respected.

For nonlinear models, sbiosteadystate uses fmincon (Optimization Toolbox) to get an initial guess for the root. The solution found by fmincon is then improved by fsolve (Optimization Toolbox). To detect convergence, sbiosteadystate uses the absolute tolerance ('AbsTol'). In other words, OptimalityTolerance, FunctionTolerance, and StepTolerance options of the corresponding optimization function are set to the 'AbsTol' value.

For linear models, sbiosteadystate finds the roots of the flux function v by solving a linear system defined by the reaction and conservation equations. For linear models, there are no rate or algebraic equations.