4. Dynamic Simulation of Models¶
This notebook example provides a basic demonstration on how to dynamically simulate models.
[1]:
import mass.test
model = mass.test.create_test_model("Phosphate_Trafficking")
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 | 0x07fe71b3f3850 |
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 = 0x7fe71b427dd0 }>
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 0x7fe71b5960b0>,
<MassSolution Phosphate_Trafficking_FluxSols at 0x7fe71b5966b0>)
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 0x7fe71b596e90>
atp_c: <scipy.interpolate.interpolate.interp1d object at 0x7fe71b596fb0>
amp_c: <scipy.interpolate.interpolate.interp1d object at 0x7fe71b5a2050>
B: <scipy.interpolate.interpolate.interp1d object at 0x7fe71b5a2110>
BP: <scipy.interpolate.interpolate.interp1d object at 0x7fe71b5a2170>
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: 0.000001
absolute_tolerance: 0.000000000001
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
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.
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 thesympy.sympify()
function. The formula can have one variable that is identical to the correspondingdict
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.40000000000001784
atp_c: 1.6000000000000718
amp_c: 0.10000000000000445
B: 1.9999999999999998
BP: 8.0
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):
[27]:
conc_sol, flux_sol = sim.find_steady_state(model, strategy="nleq2")
for metabolite, solution in conc_sol.items():
print("{0}: {1}".format(metabolite, solution))
adp_c: 0.39999999999999997
atp_c: 1.6
amp_c: 0.1
B: 2.0000000000000004
BP: 8.000000000000004
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="nleq2",
update_values=True)
model.initial_conditions # Same object as reference model in Simulation
[28]:
{<MassMetabolite adp_c at 0x7fe71b3f3950>: 0.39999999999999997,
<MassMetabolite atp_c at 0x7fe71b3f38d0>: 1.6,
<MassMetabolite amp_c at 0x7fe71b3f3c50>: 0.1,
<MassMetabolite B at 0x7fe71b3f3c10>: 2.0000000000000004,
<MassMetabolite BP at 0x7fe71b3f3c90>: 8.000000000000004}
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.2666666666666707
atp_c: 0.711111111111122
amp_c: 0.10000000000000149
B: 2.72727272727273
BP: 7.272727272727282
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: nleq2
settings:
allow_presimulation: false
presimulation_maximum_steps: 100
presimulation_time: 100
allow_approx: true
approx_tolerance: 0.000001
approx_maximum_steps: 10000
approx_time: 10000
relative_tolerance: 0.000000000001
maximum_iterations: 100
minimum_damping: 1e-20
broyden_method: 0
linearity: 3
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("maximum_iterations"))
The maximum number of iterations the solver is allowed to use (int)
For example, to change the solver options to allow for a larger maximum number of iterations:
[32]:
sim.steady_state_solver.maximum_iterations = 500
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.test.create_test_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 theSimulation
.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 0x7fe71d1201d0>,
<MassSolution Phosphate_Trafficking_Modified_ConcSols at 0x7fe71d120530>]
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 0x7fe71d120530>
For additional information on simulating multiple models, see Simulating an Ensemble of Models.