# 4. Dynamic Simulation of Models

This notebook example provides a basic demonstration on how to dynamically simulate models.

```
[1]:
```

```
import mass.example_data
model = mass.example_data.create_example_model("Phosphate_Trafficking")
```

```
Set parameter Username
```

## 4.1. Creating a Simulation

The `Simulation`

object manages all aspects related to simulating one or more model. This includes interfacing with the `RoadRunner`

object from the libRoadRunner package, which is utilized for JIT compilation of the model, integration of model ODEs, and returning simulation output.

```
[2]:
```

```
from mass import Simulation
```

### 4.1.1. Loading a model

A `Simulation`

object is initialized by providing a `MassModel`

to be loaded into the underlying `RoadRunner`

object. The provided `MassModel`

is treated as the `reference_model`

of the `Simulation`

.

`RoadRunner`

is designed to simulate models in SBML format. Therefore, models must be SBML compliant to be loaded into the `RoadRunner`

instance.

```
[3]:
```

```
sim = Simulation(reference_model=model, verbose=True)
```

```
Successfully loaded MassModel 'Phosphate_Trafficking' into RoadRunner.
```

The reference `MassModel`

is accessed using the `reference_model`

attribute.

```
[4]:
```

```
sim.reference_model
```

```
[4]:
```

Name | Phosphate_Trafficking |

Memory address | 0x07fa440bfa550 |

Stoichiometric Matrix |
5x6 |

Matrix Rank |
4 |

Number of metabolites |
5 |

Initial conditions defined |
5/5 |

Number of reactions |
6 |

Number of genes |
0 |

Number of enzyme modules |
0 |

Number of groups |
0 |

Objective expression |
0 |

Compartments |
compartment |

The underlying `RoadRunner`

instance can be accessed using the `roadrunner`

attribute:

```
[5]:
```

```
sim.roadrunner
```

```
[5]:
```

```
<roadrunner.RoadRunner() { this = 0x7fa3fa62e830 }>
```

Upon loading a model into a `Simulation`

, the numerical values for species’ initial conditions and reaction parameters are extracted. Dictionaries containing values are retrieved using the `get_model_simulation_values()`

method.

```
[6]:
```

```
initial_conditions, parameters = sim.get_model_simulation_values(model)
for metabolite, initial_condition in initial_conditions.items():
print("{0}: {1}".format(metabolite, initial_condition))
```

```
adp_c: 0.2
atp_c: 1.7
amp_c: 0.2
B: 2
BP: 8
```

## 4.2. Running Dynamic Simulations

### 4.2.1. Simulating a model

Once a model has been loaded, it can be simulated using the `Simulation.simulate()`

method. The `simulate()`

method requires a model identifier or `MassModel`

object of a loaded model and a tuple that contains the initial and final time points.

```
[7]:
```

```
# Simulate the model from 0 to 100 time units
solutions = sim.simulate(model, time=(0, 100))
solutions
```

```
[7]:
```

```
(<MassSolution Phosphate_Trafficking_ConcSols at 0x7fa460716a40>,
<MassSolution Phosphate_Trafficking_FluxSols at 0x7fa410ad4f90>)
```

After a model has been simulated, the concentration and flux solutions are returned in two specialized dictionaries known as `MassSolution`

objects.

```
[8]:
```

```
conc_sol, flux_sol = solutions
# List the first 5 points of concentration solutions
for metabolite, solution in conc_sol.items():
print("{0}: {1}".format(metabolite, solution[:5]))
```

```
adp_c: [0.2 0.20020374 0.20040726 0.20095514 0.2015016 ]
atp_c: [1.7 1.69982168 1.69964357 1.69916416 1.69868608]
amp_c: [0.2 0.19997458 0.19994916 0.19988069 0.1998123 ]
B: [2. 1.99984758 1.99969535 1.99928569 1.99887727]
BP: [8. 8.00015242 8.00030465 8.00071431 8.00112273]
```

By default, solutions are returned as `numpy.ndarrays`

. To return interpolating functions instead, the `interpolate`

argument is set as `True`

.

```
[9]:
```

```
conc_sol, flux_sol = sim.simulate(model, time=(0, 100), interpolate=True)
for metabolite, solution in conc_sol.items():
print("{0}: {1}".format(metabolite, solution))
```

```
adp_c: <scipy.interpolate._interpolate.interp1d object at 0x7fa410ae60e0>
atp_c: <scipy.interpolate._interpolate.interp1d object at 0x7fa46072c5e0>
amp_c: <scipy.interpolate._interpolate.interp1d object at 0x7fa46072c4f0>
B: <scipy.interpolate._interpolate.interp1d object at 0x7fa410ae2f40>
BP: <scipy.interpolate._interpolate.interp1d object at 0x7fa410ae6180>
```

### 4.2.2. Setting integration options

Although the integrator options are set to accomadate a variety of models, there are circumestances in which the integrator options need to be changed. The underlying integrator can be accessed using the `integrator`

property:

```
[10]:
```

```
print(sim.integrator)
```

```
< roadrunner.Integrator() >
name: cvode
settings:
relative_tolerance: 1e-06
absolute_tolerance: 1e-12
stiff: true
maximum_bdf_order: 5
maximum_adams_order: 12
maximum_num_steps: 20000
maximum_time_step: 0
minimum_time_step: 0
initial_time_step: 0
multiple_steps: false
variable_step_size: true
max_output_rows: 100000
```

Each setting comes with a brief description that can be viewed:

```
[11]:
```

```
print(sim.integrator.getDescription("variable_step_size"))
```

```
(bool) Enabling this setting will allow the integrator to adapt the size of each time step. This will result in a non-uniform time column. The number of steps or points will be ignored, and the max number of output rows will be used instead.
```

For example, to change the integration options so a uniform time vector is returned instead of one with a variable step size:

```
[12]:
```

```
sim.integrator.variable_step_size = False
# Simulate the model from 0 to 100 time units
# with output returned at evenly spaced time points
conc_sol, flux_sol = sim.simulate(model, time=(0, 100))
# Print the time vector
print(conc_sol.time[:10])
```

```
[ 0. 2. 4. 6. 8. 10. 12. 14. 16. 18.]
```

When encountering exceptions from the integrator, libRoadRunner recommends specifying an initial time step and tighter absolute and relative tolerances.

See the libRoadRunner documentation about the `roadrunner.Integrator`

class for more information on the integrator.

### 4.2.3. Simulation results and the MassSolution object

For every model simulated, two `MassSolution`

objects are returned per model. `MassSolution`

objects are always outputted as pairs, with one `MassSolution`

object containing the solutions for metabolite concentrations, and the other containing the solutions for reaction fluxes.

```
[13]:
```

```
# Simulate the model from 0 to 100 time units
sim = Simulation(model)
conc_sol, flux_sol = sim.simulate(model, time=(0, 100))
```

A `MassSolution`

for a successful simulation contains string identifiers of objects and their corresponding solutions. Because `MassSolution`

objects are specialized dictionaries, solutions can be retrieved using the object identifier as `dict`

keys. For example, to access the solution for “atp_c”:

```
[14]:
```

```
# Print first 10 solution values for ATP
conc_sol["atp_c"][:10]
```

```
[14]:
```

```
array([1.7 , 1.69982168, 1.69964357, 1.69916416, 1.69868608,
1.69820944, 1.69773427, 1.69673394, 1.69460788, 1.6925119 ])
```

If care is taken when assigning object identifiers (e.g., does not start with a number, does not contain certain characters such as “-”), it is possible to access solutions inside of a `MassSolution`

as if the corresponding keys were attributes.

```
[15]:
```

```
# Print first 10 solution values for ATP
conc_sol.atp_c[:10]
```

```
[15]:
```

```
array([1.7 , 1.69982168, 1.69964357, 1.69916416, 1.69868608,
1.69820944, 1.69773427, 1.69673394, 1.69460788, 1.6925119 ])
```

The time points returned by the integrator are accessible using the `MassSolution.time`

attribute:

```
[16]:
```

```
# Print the first 10 time points
print(conc_sol.time[0:10])
```

```
[0.00000000e+00 8.47847116e-08 1.69569423e-07 3.98257356e-07
6.26945288e-07 8.55633221e-07 1.08432115e-06 1.56811392e-06
2.60709564e-06 3.64607737e-06]
```

The solutions contained within the `MassSolution`

can be obtained as a `pandas.DataFrame`

using the `to_frame()`

method.

```
[17]:
```

```
conc_sol.to_frame()
```

```
[17]:
```

adp_c | atp_c | amp_c | B | BP | |
---|---|---|---|---|---|

Time | |||||

0.000000e+00 | 0.200000 | 1.700000 | 0.200000 | 2.000000 | 8.000000 |

8.478471e-08 | 0.200204 | 1.699822 | 0.199975 | 1.999848 | 8.000152 |

1.695694e-07 | 0.200407 | 1.699644 | 0.199949 | 1.999695 | 8.000305 |

3.982574e-07 | 0.200955 | 1.699164 | 0.199881 | 1.999286 | 8.000714 |

6.269453e-07 | 0.201502 | 1.698686 | 0.199812 | 1.998877 | 8.001123 |

... | ... | ... | ... | ... | ... |

3.012599e+01 | 0.399998 | 1.599991 | 0.099999 | 2.000000 | 8.000000 |

3.828512e+01 | 0.399998 | 1.599992 | 0.100000 | 2.000000 | 8.000000 |

5.517474e+01 | 0.399998 | 1.599994 | 0.100000 | 2.000000 | 8.000000 |

8.169531e+01 | 0.399999 | 1.599996 | 0.100000 | 2.000000 | 8.000000 |

1.000000e+02 | 0.399999 | 1.599997 | 0.100000 | 2.000000 | 8.000000 |

171 rows × 5 columns

Solutions also can be viewed visually using the `view_time_profile()`

method. Note that this requires `matplotlib`

to be installed in the environment. See Plotting and Visualization for more information.

```
[18]:
```

```
conc_sol.view_time_profile()
```

### 4.2.4. Aggregate variables and solutions

Often, it is desirable to look at mathematical combinations of metabolites concentrations or reaction fluxes. To create an aggregate variable, the `MassSolution.make_aggregate_solution()`

method is used. To use the method, three inputs are required:

A unique ID for the aggregate variable.

The mathematical equation for the aggregate variable given as a

`str`

.A list of the

`MassSolution`

keys representing variables used in the equation.

For example, to make the Adenylate Energy Charge [Atk68], the occupancy and capacity pools are first defined:

```
[19]:
```

```
occupancy = conc_sol.make_aggregate_solution(
aggregate_id="occupancy",
equation="(atp_c + 0.5 * adp_c)",
variables=["atp_c", "adp_c"])
conc_sol.update(occupancy)
print(list(conc_sol.keys()))
```

```
['adp_c', 'atp_c', 'amp_c', 'B', 'BP', 'occupancy']
```

The aggregate variables are returned as a `dict`

, which can be added to the `MassSolution`

object. Alternatively, the `update`

flag can be set as `True`

to automatically add an aggregate variable to the solution after creation.

```
[20]:
```

```
capacity = conc_sol.make_aggregate_solution(
aggregate_id="capacity",
equation="(atp_c + adp_c + amp_c)",
variables=["atp_c", "adp_c", "amp_c"],
update=True)
print(list(conc_sol.keys()))
```

```
['adp_c', 'atp_c', 'amp_c', 'B', 'BP', 'occupancy', 'capacity']
```

Aggregate variables formed from other aggregate variables also can be created using the `make_aggregate_solution()`

method as long as the aggregate variables have been added to the `MassSolution`

. To make the energy charge from the occupancy and capacity aggregate variables:

```
[21]:
```

```
ec = conc_sol.make_aggregate_solution(
aggregate_id="energy_charge",
equation="occupancy / capacity",
variables=["occupancy", "capacity"],
update=True)
```

If care is taken when assigning aggregate variable identifiers, it is possible to access aggregate variable solutions inside of a `MassSolution`

, as if aggregate variable keys were attributes.

```
[22]:
```

```
# Print first 10 solution points for the energy charge
conc_sol.energy_charge[:10]
```

```
[22]:
```

```
array([0.85714286, 0.85710645, 0.8570701 , 0.85697226, 0.85687471,
0.85677749, 0.85668059, 0.85647669, 0.85604374, 0.85561747])
```

## 4.3. Perturbing a Model

To simulate various disturbances in the system, the `perturbations`

argument of the `simulate()`

method can be used. There are several types of perturbations that can be implemented for a given simulation as long as they adhere to the following guidelines:

Perturbations are provided to the method as a

`dict`

with dictionary keys that correspond to variables to be changed. Dictionary values are the new numerical values or mathematical expressions as strings that indicate how the value is to be changed.A formula for the perturbation can be provided as a

`str`

as long as the formula string can be sympified via the`sympy.sympify()`

function. The formula can have one variable that is identical to the corresponding`dict`

key.Boundary conditions can be set as a function of time. The above rules still apply, but allow for the time “t”, as a second variable.

Some examples are demonstrated below.

A simulation without perturbations:

```
[23]:
```

```
conc_sol, flux_sol = sim.simulate(model, time=(0, 1000))
conc_sol.view_time_profile()
```

Perturbing the initial concentration of ATP from 1.6 to 2.5:

```
[24]:
```

```
conc_sol, flux_sol = sim.simulate(
model, time=(0, 1000), perturbations={"atp_c": 2.5})
conc_sol.view_time_profile()
```

Increasing the rate constant of ATP use by 50%:

```
[25]:
```

```
conc_sol, flux_sol = sim.simulate(
model, time=(0, 1000), perturbations={"kf_use": "kf_use * 1.5"})
conc_sol.view_time_profile()
```

## 4.4. Determining Steady State

The steady state for models can be found using the `Simulation.find_steady_state()`

method. This method requires a model identifier or a `MassModel`

object and a string thats indicates a strategy for finding the steady state. For example, to find the steady state by simulating the model for a long time:

```
[26]:
```

```
sim = Simulation(reference_model=model)
conc_sol, flux_sol = sim.find_steady_state(model, strategy="simulate")
for metabolite, solution in conc_sol.items():
print("{0}: {1}".format(metabolite, solution))
```

```
adp_c: 0.3999999999999983
atp_c: 1.5999999999999932
amp_c: 0.09999999999999959
B: 1.999999999999997
BP: 7.999999999999988
```

Alternatively, a steady state solver can be utilized with a root-finding algorithm to determine the steady state. For example, to use a non-linear equation solver that implements a global Newton method with adaptive damping strategies (NLEQ2) or a line-search method:

```
[27]:
```

```
conc_sol, flux_sol = sim.find_steady_state(model, strategy="newton_linesearch")
for metabolite, solution in conc_sol.items():
print("{0}: {1}".format(metabolite, solution))
```

```
adp_c: 0.39999920693137425
atp_c: 1.599996808967697
amp_c: 0.09999980294303153
B: 2.0000000187666185
BP: 7.999999981233392
```

Setting `update_values=True`

updates the model initial conditions and fluxes with the steady state solution:

```
[28]:
```

```
conc_sol, flux_sol = sim.find_steady_state(model, strategy="newton_linesearch",
update_values=True)
model.initial_conditions # Same object as reference model in Simulation
```

```
[28]:
```

```
{<MassMetabolite adp_c at 0x7fa4606bf730>: 0.39999920693137425,
<MassMetabolite atp_c at 0x7fa4606bf6a0>: 1.599996808967697,
<MassMetabolite amp_c at 0x7fa4606bf700>: 0.09999980294303153,
<MassMetabolite B at 0x7fa4606bf790>: 2.0000000187666185,
<MassMetabolite BP at 0x7fa4606bf7c0>: 7.999999981233392}
```

The `find_steady_state()`

method also allows for perturbations to be made before determining a steady state solution:

```
[29]:
```

```
conc_sol, flux_sol = sim.find_steady_state(
model, strategy="simulate", perturbations={"kf_use": "kf_use * 1.5"})
for metabolite, solution in conc_sol.items():
print("{0}: {1}".format(metabolite, solution))
```

```
adp_c: 0.26666666666666217
atp_c: 0.7111111111110989
amp_c: 0.09999999999999833
B: 2.72727272727273
BP: 7.272727272727279
```

### 4.4.1. The steady state solver

Although the steady state solver options are set to accommodate a variety of models, there are circumstances in which the steady state solver options need to be changed. The underlying solver can be accessed using the `steady_state_solver`

property:

```
[30]:
```

```
print(sim.steady_state_solver)
```

```
< roadrunner.SteadyStateSolver() >
name: newton_linesearch
settings:
auto_moiety_analysis: false
allow_presimulation: true
presimulation_time: 5
presimulation_times: [0.1, 1, 10, 100, 1000, 10000]
presimulation_maximum_steps: 500
allow_approx: true
approx_tolerance: 1e-06
approx_maximum_steps: 10000
approx_time: 10000
num_max_iters: 200
allow_negative: false
print_level: 0
eta_form: 'eta_choice1'
no_init_setup: false
no_res_monitoring: false
max_setup_calls: 10
max_subsetup_calls: 5
eta_constant_value: 0.1
eta_param_gamma: 0
eta_param_alpha: 0
res_mon_min: 1e-05
res_mon_max: 0.9
res_mon_constant_value: 0.9
no_min_eps: false
max_newton_step: 0
max_beta_fails: 10
func_norm_tol: 0
scaled_step_tol: 0
rel_err_func: 0
```

Analogous to the integrator, each setting for the steady state solver comes with a brief description that can be viewed:

```
[31]:
```

```
print(sim.steady_state_solver.getHint("approx_maximum_steps"))
```

```
Maximum number of steps that can be taken for steady state approximation routine (int).
```

For example, to change the solver options to allow for a larger maximum number of iterations:

```
[32]:
```

```
sim.steady_state_solver.approx_maximum_steps = 15000
```

For more information on steady state solver options, see the libRoadRunner documentation about the `roadrunner.SteadyStateSolver`

class.

## 4.5. Simulating Multiple Models

Multiple models can be added to a `Simulation`

object in order to perform simulations on several models. Below, a simple example utilizing a copy of the model with a smaller total buffer concentration is made for demonstration purposes:

```
[33]:
```

```
model = mass.example_data.create_example_model("Phosphate_Trafficking")
sim = Simulation(model)
# Make a modified model
modified = model.copy()
modified.id = "Phosphate_Trafficking_Modified"
modified.update_initial_conditions({"BP": 4, "B": 1})
```

To add an additional model to an existing `Simulation`

object, three criteria must be met:

The model must have equivalent ODEs to the

`reference_model`

used in creating the`Simulation`

.All models in the

`Simulation`

must have unique identifiers.Numerical values necessary for simulation must be already defined for a model.

Use the `MassModel.has_equivalent_odes()`

method to check if ODEs are equivalent.

```
[34]:
```

```
sim.reference_model.has_equivalent_odes(modified, verbose=True)
```

```
[34]:
```

```
True
```

Use the `Simulation.add_models()`

method to load the additional model.

```
[35]:
```

```
sim.add_models(models=[modified], verbose=True)
```

```
Successfully loaded MassModel 'Phosphate_Trafficking_Modified'.
```

Use the `simulate()`

method to simulate multiple models by providing a list of model objects or their identifiers.

```
[36]:
```

```
conc_sol_list, flux_sol_list = sim.simulate(
models=[model, modified], time=(0, 100))
```

After simulating multiple models, the `MassSolution`

objects are returned in two `cobra.DictList`

objects. The first `DictList`

contains the concentrations solutions, and the second `DictList`

contains the flux solutions for simulated models.

```
[37]:
```

```
conc_sol_list
```

```
[37]:
```

```
[<MassSolution Phosphate_Trafficking_ConcSols at 0x7fa400e647c0>,
<MassSolution Phosphate_Trafficking_Modified_ConcSols at 0x7fa46072c360>]
```

The `get_by_id()`

method can be used to access a specific solution:

```
[38]:
```

```
conc_sol_list.get_by_id("_".join((modified.id, "ConcSols")))
```

```
[38]:
```

```
<MassSolution Phosphate_Trafficking_Modified_ConcSols at 0x7fa46072c360>
```

For additional information on simulating multiple models, see Simulating an Ensemble of Models.