Documentation

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

example

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

example

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

example

``````[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, `[]`.```

example

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

example

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

example

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

## Examples

collapse all

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

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

collapse all

SimBiology model, specified as a SimBiology `Model 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.

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 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'`

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(‖\frac{d\stackrel{\to }{S}}{dt}‖, where $\stackrel{\to }{S}$ is a vector of nonconstant species, parameters, and compartments.

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(‖\frac{d\stackrel{\to }{S}}{dt}‖, where $\stackrel{\to }{S}$ is a vector of non-constant species, parameters, and compartments.

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.

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

collapse all

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

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.

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.

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.

collapse all

### 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` to get an initial guess for the root. The solution found by `fmincon` is then improved by `fsolve`. 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.