3. Dynamic Simulation: The Basic Procedure

Once a set of dynamic mass balance equations has been formulated, they can be numerically solved, and thus the behavior of a network can be simulated in response to environmental and genetic changes. Simulation results can be obtained using a number of different software packages. Dynamic simulation generates the time dependent behavior of the concentrations, i.e., \(\textbf{x}\)(t). This solution can be obtained in response to several different types of perturbations and the results graphically displayed. The basic principles and procedures associated with dynamic simulation are covered in this chapter. The following three chapters then apply the simulation process to a set of simple but progressively more complex and relevant examples.

3.1. Numerical Solutions

Network dynamics are described by a set of ordinary differential equations (ODEs): the dynamic mass balance equations; see Eq. (1.1). To obtain the dynamic solutions, we need three things: first, the equations themselves; second, the numerical values for the kinetic constants that are in the equations; and third, the initial conditions and parameters that are being perturbed. We describe each briefly.

1. To formulate the mass balances we have to specify the system boundary, the fluxes in and out of the system, and the reactions that take place in the network. From the set of reactions that are taking place, a stoichiometric matrix is formed. This matrix is then put into Eq. (1.1) . One can multiply out the individual dynamic mass balances, as was done in Eq. (2.13) for the adenosine phosphate network, to prevent a large number of numerical operations that involve multiplication of reaction rates by the zero elements in \(\textbf{S}\). The reaction rate laws for the reactions are then identified and substituted into the equations. Typically, one would use elementary kinetics as shown in Eq. (2.6), or apply more complex rate laws if they are appropriate and available. This process leads to the definition of the dynamic mass balances.

2. The numerical values for the kinetic parameters in the rate laws have to be specified, as do any imposed fluxes across the system boundary. Obtaining numerical values for the kinetic constants is typically difficult. They are put into a parameter vector designated by \(\textbf{k}\). In select cases, detailed kinetic characterization has been carried out. More often, though, one only knows these values approximately. It is important to make sure that all units are consistent throughout the equations and that the numerical values used are in the appropriate units.

3. With the equations and numerical values for the kinetic constants specified \((\textbf{k})\), we can simulate the responses of the network that they represent. To do so, we have to set initial conditions \((x_0)\). This leads to the numerical solution of

\[\begin{equation} \frac{d\textbf{x}}{dt} = \textbf{Sv(x;k)},\ \textbf{x}(t = 0) = \textbf{x}_0 \tag{3.1} \end{equation}\]

There are three conditions that are typically considered.

1. First, the initial conditions for the concentrations are set, and the motion of the network into a steady state (open system) or equilibrium (closed system) is simulated. This scenario is typically physiologically unrealistic since individual concentrations cannot simply change individually in a living cell.

2. Second, a change in an input flux is imposed on a network that is in a steady state. This scenario can be used to simulate the response of a cell to a change in its environment.

3. Third, a change in a kinetic parameter is implemented at the initial time. The initial concentrations are typically set at the steady state values with the nominal value of the parameter. The equations are then simulated to a long time to obtain the steady state values that correspond to the altered kinetic parameters. These are set as the initial conditions when examining the responses of the system with the altered kinetic properties.

4. Once the solution has been obtained it can be graphically displayed and the results analyzed. There are several ways to accomplish this step, as detailed in the next two sections. The analysis of the results can lead to post-processing of the output to form an alternative set of dynamic variables.

The simulation is implemented using a numerical solver. Currently, such implementation is carried out using standard and readily available software, such as Mathematica or MATLAB. Specialized simulation packages are also available (Table 3.1). After the simulation is set up and the conditions specified, the software computes the concentrations as a function of time. The output is a file that contains the numerical values of the concentrations at a series of time points (Figure 3.1). This set of numbers is typically graphically displayed, and/or used for subsequent computations.

Table 3.1: Available software for dynamic simulation. Assembled by Neema Jamshidi.

We will be using iPython notebooks using the python software package MASSpy as our simulation software.

3.2. Graphically Displaying the Solution

The simulation procedure described in the previous section results in a file that contains the concentrations as a function of time (Figure 3.1). These results are graphically displayed, typically in two ways: by plotting the concentrations as a function of time, or by plotting two concentrations against one another with time as a parameter along the trajectory.

Figure 3.1: The fundamental structure of the file \(\textbf{x}\)(t) that results from a numerical simulation. The two vertical bars show the list of values that would be used to compute \(\sigma_{12}\)(2) (see Eq. 3.8); that is, the correlation between \(x_1\) and \(x_2\) with a time lag of 2.

Before describing these methods, we observe certain fundamental aspects of the equations that we are solving. The dynamic mass balances can be expanded as:

\[\begin{equation} \frac{d\textbf{x}}{dt} = \textbf{Sv(x)} = \sum\textbf{s}_i v_i(\textbf{x}) \tag{3.2} \end{equation}\]

In other words, the time derivatives are linear combinations of the reaction vectors \((\textbf{s}_i)\) weighted by the reaction rates, that in turn change with the concentrations that are time varying. Thus, the motions are linear combinations of the directions specified by \(\textbf{s}_i\). This characteristic is important because if the \(v_i\). have different time constants, the motion can be decomposed in time along these reaction vectors.

3.2.1. Time Profiles

The simulation results in a file that contains the vector \(\textbf{x}\)(t) and the time points at which the numerical values for the concentrations are given. These time points can be specified by the user or are automatically generated by the solver used. Typically, the user specifies the initial time, the final time, and sometimes the time increment between the time points where the simulator stores the computed concentration values in the file. The results can then be graphically displayed depending on a few features of the solution. Some of these are shown in Figure 3.2 and are now briefly described:

• Panel A: The most common way to display a dynamic solution is to plot the concentration as a function of time.

• Panel B: If there are many concentration variables they are often displayed on the same graph.

• Panel C: In many cases there are different response times and one plots multiple time profiles where the x-axis on each plot is scaled to a particular response time. Alternatively, one can use a logarithmic scale for time.

• Panel D: If a variable moves on many time scales changing over many orders of magnitude, the y-axis is often displayed on a logarithmic scale.

Figure 3.2: Graphing concentrations over time. (a) A single concentration shown as a function of time. (b) Many concentrations shown as a function of time. (c) A single concentration shown as a function of time separately on different time scales. (d) The logarithm of a single concentration shown as a function of time to distinguish the decay on different time scales.

The solution can thus be displayed in different ways depending on the characteristics of the time profiles. One normally plays with these representations to get an understanding of the responses of the network that they have formulated and to represent the features in which one is interested.

3.2.2. Dynamic phase portraits

Dynamic phase portraits represent trajectories formed when two concentrations plotted against each other, parameterized with respect to time (Figure 3.3). The dynamic trajectories in the diagram move from an initial state to a final state. Analysis of these trajectories can point to key dynamic relationships between compounds in a biochemical reaction network. For example, if a system is dynamically stable, the dynamic trajectories will converge to a single point in the plane, known as an attracting fixed pointattracting fixed point. A stable steady-state point would represent a homeostatic state. Conversely, if the system is unstable, the trajectories will not approach a fixed point but diverge away from it. The former is essentially always the case for biochemical reaction networks representing real cells. The way the trajectories converge on the steady state is highly informative as different dynamic characteristics are evident from the trajectory.

Figure 3.3: A dynamic phase portrait.

3.2.3. Characteristic features of phase portraits

A trajectory in the phase portraitphase portrait may indicate the presence of one or more general dynamic features. Namely, the shapes of the trajectories contain significant information about the dynamic characteristics of a network. Some important features of trajectories in a phase portrait are shown in Figure 3.4

Figure 3.4: General features of dynamic phase portraits. Dynamic phase portraits are formed by graphing the time dependent concentrations of two concentrations \((x_1\) and \(x_2)\) against one another. Phase portraits have certain characteristic features. (a) Conservation relationship. (b) A pair of concentrations that could be in quasi-equilibrium with one another. (c) Motion of the two concentrations dynamically independent of one another. (d) Closed loop traces representing either a periodic motion or a return to the original steady state. Modified from Kauffman 2002 [64].

1. When the trajectory has a negative slope, it indicates that one concentration is increasing while the other is decreasing. The concentrations are moving on the same time scales but in opposite directions; that is, one is consumed while the other is produced. This feature might represent the substrate concentration versus the product concentration of a given reaction. Such behavior helps define aggregate concentration variablesaggregate concentration variables.

2. When a trajectory in the phase portrait between two concentrations is a straight line with a positive slope, it means that the two concentrations are moving in tandem; i.e., as one increases so does the other. This feature is observed when two or more concentrations move on the same time scales and are in quasi-equilibrium with one another. Such behavior helps define aggregate concentration variables.

3. When a trajectory is vertical or horizontal, it indicates that one of the concentrations is changing while the other remains constant. This feature implies either that the motions of the concentrations during the trajectory are independent of one another or that the dynamic motions of the concentrations progress on different characteristic time scales. Such behavior helps define time scale decomposition.

4. When a trajectory forms a closed loop, it implies one of two possibilities. The system never converges to a steady state over time but oscillates forming a closed loop trajectory. On the other hand, if the orbit begins at one point, moves away from it, then returns to the same point after a sufficiently long time interval, then it implies that a change in another variable in the system forced it away from its steady state temporarily, but it returned to the original steady state. Such behavior helps define disturbance rejection characteristics.

Figure 3.5: A schematic of a tiled phase portrait.The matrix is symmetric, making it possible to display statistical information about a phase portrait in the mirror position.The diagonal elements are meaningless.Originally developed in Kauffman 2002 [64].

The qualitative characteristics of dynamic phase portraitsphase portrait can provide insight into the dynamic features of a network. A trajectory may have more than one of these basic features. For instance, there can be a fast independent motion (i.e., a horizontal phase portrait trajectory) followed by a line with a positive slope after an equilibrium state has been reached.

3.2.4. Tiling dynamic phase portraits

Phase portraits show the dynamic relationships between two variables on multiple time scales, see Figure 3.5. If a system has \(n\) variables, then there are \(n^2\) dynamic phase portraits. All pair-wise phase portraits can be tiled in a matrix form where the \(\textit{i}\), \(\textit{j}\) entry represents the dynamic phase portrait between variables \(x_i\) and \(x_j\). Note that such an array is symmetric and that the diagonal elements are un-informative. Thus, there are \((n^2-n)/2\) phase portraits of interest. This feature of this graphical representation opens the possibility of putting the phase portrait in the \(\textit{i}\), \(\textit{j}\) position in the array and showing other information (such as a regression coefficient or a slope) in the corresponding \(\textit{j}\), \(\textit{i}\) position.

Since the time scales in biochemical reaction networks typically vary over many orders of magnitude, it often makes sense to make a series of tiled phase portraits, each of which represents a key time scale. For instance, rapid equilibration leads to straight lines with positive slopes in the phase portrait (Figure 3.4b) where the slope is the equilibrium constant of the reaction. This may be one of many dynamic events taking place. If a phase portrait is graphed separately on this time scale alone, the positive line will show up with a high regression coefficient and a slope that corresponds to the equilibrium constant.

3.3. Post-Processing the Solution

The initial suggestions obtained from graphing and visualizing the concentration vector \(\textbf{x}\)(t) can lead to a more formal analysis of the results. We describe three post-processing procedures of \(\textbf{x}\)(t).

3.3.1. Computing the fluxes from the concentration variables:

The solution for the concentrations \(\textbf{x}\)(t)can be used to compute the fluxes from

\[\begin{equation} \textbf{v}(t)= \textbf{v}(\textbf{x}(t)) \tag{3.3} \end{equation}\]

and subsequently we can plot the fluxes in the same way as the concentrations. Graphical information about both the \(\textbf{x}\)(t) and \(\textbf{v}\)(t) is useful.

3.3.2. Combining concentrations to form aggregate variables:

The graphical and statistical multi-time scale analysis discussed above may lead to the identification of aggregate variables. Pooled variables, p, are computed from

\[\begin{equation} \textbf{p}(t)= \textbf{Px}(t)) \tag{3.4} \end{equation}\]

where the pool transformation matrix, \(\textbf{P}\), defines the linear combination of the concentration variables that forms the aggregate variables. For instance, if we find that a logical way to pool two variables, \(x_1\) and \(x_2\), into new aggregate variables is \(p_1 = x_1 + x_2\) and \(p_2 = x_1 - x_2\), then we form the following matrix equation describing these relationships as:

\[\begin{split}\begin{equation} \textbf{p}(t) = \textbf{Px}(t) = \begin{pmatrix} {p_1(t)} \\ {p_2(t)} \end{pmatrix} = \begin{pmatrix} {1} & {1} \\ {1} & {-1} \end{pmatrix} \begin{pmatrix} {x_1(t)} \\ {x_2(t)} \end{pmatrix} = \begin{pmatrix} {x_1(t) + x_2(t)} \\ {x_1(t) - x_2(t)} \end{pmatrix} \end{equation}\end{split}\]

The dynamic variables, \(\textbf{p}\)(t), can be graphically studied as described in the previous section.

3.3.2.1. Example: The Phosphorylated Adenosines

The pool formation discussed in Chapter 2 can be described by the pool transformation matrix:

\[\begin{split}\begin{equation}\textbf{P} = \begin{pmatrix} {1} & {1} & {0} \\ {2} & {1} & {0} \\ {1} & {1} & {1} \end{pmatrix} \end{equation}\end{split}\]

\[\tag{3.5}\]

and thus

\[\begin{split}\begin{equation}\\textbf{p} = \textbf{Px} = \textbf{P}\begin{pmatrix} {\text{ATP}} \\ {\text{ADP}} \\ {\text{AMP}} \end{pmatrix} = \begin{pmatrix} {\text{ATP} + \text{ADP}} \\ {2 \text{ATP} + \text{ADP}} \\ {\text{ATP} + \text{ADP} + \text{AMP}} \end{pmatrix}\end{equation}\end{split}\]

\[\tag{3.6}\]

The pool sizes \(p_i\)(t) can then be graphed as a function of time.

3.3.3. Correlating concentrations over time:

One can construct the time-separated correlation matrix, \(\textbf{R}\), based on a time scale structure of a system. In this matrix, we compute the correlation between two concentrations on a time scale as:

\[\begin{equation} \textbf{R}(\tau) = (r_{ij}) = \frac{\sigma_{ij}(\tau)}{\sqrt{\sigma_{ii}\sigma_{jj}}} \tag{3.7} \end{equation}\]

in which \(\sigma_{ii}\) is the variance of the dataset \(x_i(k)\) and \(\sigma_{ij}(\tau)\) is the time-lagged covariance between the discrete, uniformly sampled datasets \(x_i(k)\) and \(x_j(k + \tau)\), determined as,

\[\begin{equation} \sigma_{ij}(\tau) = \frac{1}{n}\sum\limits_{k=1}^{n-\tau} (x_i(k) - \bar{x_i})(x_j(k + \tau) - \bar{x_j}) \tag{3.8} \end{equation}\]

in which \(n\) is the number of data points in the series, and \(\bar{x_i}\) indicates the average value of the series \(x_i\). The values in \(\textbf{R}\) range from -1 to 1, indicating perfect anti-correlation or correlation, respectively, between two datasets with a delay of time steps. Elements in \(\textbf{R}\) equal to zero indicate that the two corresponding datasets are completely uncorrelated. If such correlation computations were done for the cases shown in Figure 3.4, one would expect to find a strong negative correlation for the data shown in Figure 3.4a, a strong positive correlation for Figure 3.4b, and no correlation for Figure 3.4c,

The correlation computations can be performed with an increment, \(\tau\), offset in time between two concentrations. An example of a time offset is shown in Figure 3.2 showing the values used from the output file to compute the correlation between \(x_1\) and \(x_2\) with a time lag of 2.

The matrix of phase portraits is symmetric with uninformative diagonal elements. One can therefore enter a correlation coefficient corresponding to a particular phase portrait in the transpose position to the phase portrait in the matrix. A correlation coefficient provides a quantitative description of the phase portrait’s linearity between the two variables over the time scale displayed. In addition to the correlation coefficient, the slope can be computed and displayed, giving the equilibrium constant between the two compounds displayed.

3.4. Demonstration of the Simulation Procedure in MASSpy

3.4.1. Setting up the model

The following builds the model of three reactions in series that is described on pages 51-56 in the book. We show how the model is built, simulated, solutions graphically displayed, solutions post processed and analyzed mathematically.

To construct a model in MASSpy, the MassModel, MassReaction, and MassMetabolite objects need to be imported into the environment.

from mass import MassModel, MassMetabolite, MassReaction

3.4.1.1. Defining metabolites and reactions

One method for creating the model is to objects that represent the metabolites and reactions. Metabolite are represented by MassMetabolite objects, and can be created by providing a unique identifier for that object. Therefore we can define the four metabolites, \(x_1, x_2, x_3\), and \(x_4\) by the following;

x1 = MassMetabolite('x1')
x2 = MassMetabolite('x2')
x3 = MassMetabolite('x3')
x4 = MassMetabolite('x4')

Reactions are represented by MassReaction objects, and like metabolites, they can be also created by providing a unique identifier for that object.

v1 = MassReaction('v1')
v2 = MassReaction('v2')

By default, a reaction is considered reversible. However, if we wish to make an irreversible reaction, we set the reversible argument to False.

v3 = MassReaction('v3', reversible=False)

Once the MassReaction objects have been created, metabolites can be added to the reaction using the MassReaction.add_metabolites method. To quickly see how this method is used, we can use the help() function. Alternatively, we can go to the API documentation and read about how the MassReaction.add_metabolites method works.

help(MassReaction.add_metabolites)

Help on function add_metabolites in module mass.core.mass_reaction:

add_metabolites(self, metabolites_to_add, combine=True, reversibly=True)
    Add metabolites and their coefficients to the reaction.

    If the final coefficient for a metabolite is 0 then it is removed from
    the reaction.

    The change is reverted upon exit when using the :class:`~.MassModel`
    as a context.

    Notes
    -----
    * A final coefficient of < 0 implies a reactant and a final
      coefficient of > 0 implies a product.

    * Extends :meth:`~cobra.core.reaction.Reaction.add_metabolites` of the
      :class:`cobra.Reaction <cobra.core.reaction.Reaction>` by first
      ensuring that the metabolites to be added are
      :class:`.MassMetabolite`\ s and not
      :class:`cobra.Metabolites <cobra.core.metabolite.Metabolite>`.
      and error message raised reflects the :mod:`mass` object.

    * If a :class:`cobra.Metabolite <cobra.core.metabolite.Metabolite>` is
      provided. a warning is raised and a :class:`.MassMetabolite`
      will be instantiated using the
      :class:`cobra.Metabolite <cobra.core.metabolite.Metabolite>`.

    Parameters
    ----------
    metabolites_to_add : dict
        A ``dict`` with :class:`.MassMetabolite`\ s or metabolite
        identifiers as keys and stoichiometric coefficients as values. If
        keys are strings (id of a metabolite), the reaction must already
        be part of a :class:`~.MassModel` and a metabolite with the given
        id must already exist in the :class:`~.MassModel`.
    combine : bool
        Describes the behavior of existing metabolites.
        If ``True``, the metabolite coefficients are combined together.
        If ``False`` the coefficients are replaced.
    reversibly : bool
        Whether to add the change to the context to make the change
        reversible (primarily intended for internal use).

    See Also
    --------
    :meth:`subtract_metabolites`

To use MassReaction.add_metabolites, a dictionary input is required, where the MassMetabolite objects are keys and the value is their stoichiometric coefficient. Reactants are defined with negative coefficients, while products are defined with positive coefficients.

v1.add_metabolites({x1 : -1, x2 : 1})
v2.add_metabolites({x2 : -1, x3 : 1})
v3.add_metabolites({x3 : -1, x4 : 1})

Reactions, e.g., \(v_1\) can be used to define any kind of chemical transformation, association, activation etc. A series of methods are provided for inspection of the reaction.

v1.id

'v1'

v1.reactants

[<MassMetabolite x1 at 0x7fd978e75190>]

v1.products

[<MassMetabolite x2 at 0x7fd978e75070>]

v1.stoichiometry

[-1, 1]

v1.reversible

True

Check the documentation for the MassReaction class for further details.

3.4.1.2. Model Setup

To construct a model capable of dynamic simulation, a MassModel object must be created. The minimal input for creating a MassModel object is a unique identifier.

model = MassModel('Model')

Set parameter Username

To add reactions and their corresponding metabolites to the model, the MassModel.add_reactions method can be used by providing a list of reactions to add to the model.

model.add_reactions([v1, v2, v3])

3.4.1.3. Model Inspection

Similar to the MassReaction object, the MassModel object also has various methods that can be used to inspect the model. For example, to obtain the list of reactions and species in the system:

model.reactions

[<MassReaction v1 at 0x7fd978e755b0>,
 <MassReaction v2 at 0x7fd978e752e0>,
 <MassReaction v3 at 0x7fd978e75760>]

model.metabolites

[<MassMetabolite x1 at 0x7fd978e75190>,
 <MassMetabolite x2 at 0x7fd978e75070>,
 <MassMetabolite x3 at 0x7fd978e75340>,
 <MassMetabolite x4 at 0x7fd978e756d0>]

In some circumstances, it is helpful to iterate through a reaction and its associated metabolites using a loop:

print("Model ID: %s" % model.id)
for rxn in model.reactions:
    print("\nReaction: %s\n------------" % rxn.id)
    for metab, stoichiometry in rxn.metabolites.items():
        print("%s: %s " % (metab.id, stoichiometry))

Model ID: Model

Reaction: v1
------------
x1: -1
x2: 1

Reaction: v2
------------
x3: 1
x2: -1

Reaction: v3
------------
x4: 1
x3: -1

To examine the stoichiometric matrix:

model.S

array([[-1.,  0.,  0.],
       [ 1., -1.,  0.],
       [ 0.,  1., -1.],
       [ 0.,  0.,  1.]])

The stoichiometric matrix can also be viewed as a pandas.DataFrame with annotated information about the metabolites and reactions.

Note: The update_model argument can be used to store matrix as the specified array_type for the next time the stoichiometric matrix is viewed.

model.update_S(array_type="DataFrame", update_model=True)

	v1	v2	v3
x1	-1.0	0.0	0.0
x2	1.0	-1.0	0.0
x3	0.0	1.0	-1.0
x4	0.0	0.0	1.0

The rate equations can be examined,

for rxn, rate in model.rates.items():
    print("%s: %s" % (rxn.id, rate))

v1: kf_v1*(x1(t) - x2(t)/Keq_v1)
v2: kf_v2*(x2(t) - x3(t)/Keq_v2)
v3: kf_v3*x3(t)

or just one rate equation can be called out:

print(model.rates[v2])

kf_v2*(x2(t) - x3(t)/Keq_v2)

The ordinary differential equations can be also be listed in full,

for metab, ode in model.odes.items():
    print("%s: %s" % (metab.id, ode))

x1: -kf_v1*(x1(t) - x2(t)/Keq_v1)
x2: kf_v1*(x1(t) - x2(t)/Keq_v1) - kf_v2*(x2(t) - x3(t)/Keq_v2)
x3: kf_v2*(x2(t) - x3(t)/Keq_v2) - kf_v3*x3(t)
x4: kf_v3*x3(t)

or just one ordiniary differential equation can be called out:

print(model.odes[x3])

kf_v2*(x2(t) - x3(t)/Keq_v2) - kf_v3*x3(t)

Note that none of these expressions have been provided during the model construction process. Instead the expresions have been generated automatically from the provided list of reactions and their metabolites.

3.4.1.4. Set parameters and initial condtions

When using Jupyter notebooks, an overview of the model is rendered as a table when only the model object is called. Note that this also applies to metabolites and reactions.

model

Name	Model
Memory address	0x07fd9893fa490
Stoichiometric Matrix	4x3
Matrix Rank	3
Number of metabolites	4
Initial conditions defined	0/4
Number of reactions	3
Number of genes	0
Number of enzyme modules	0
Number of groups	0
Objective expression	0
Compartments

From the model overview it can be seen that no parameters or initial conditions have been defined. Parameters can be defined directly for a specific reaction:

v1.forward_rate_constant = 1
v2.kf = 0.01 # Shorthand method
v3.kf = 0.0001

v1.equilibrium_constant = 1
v2.Keq = 1 # Shorthand method

for param_type, param_dict in model.parameters.items():
    print("%s: %s" %(param_type, param_dict))

kf: {'kf_v1': 1, 'kf_v2': 0.01, 'kf_v3': 0.0001}
Keq: {'Keq_v1': 1, 'Keq_v2': 1, 'Keq_v3': inf}
kr: {}
v: {}
Custom: {}
Boundary: {}

Initial conditions for metabolites can be defined directly for a specific metabolite,

x1.initial_condition = 1
x2.ic  = 0 # Shorthand method
model.initial_conditions

{<MassMetabolite x1 at 0x7fd978e75190>: 1,
 <MassMetabolite x2 at 0x7fd978e75070>: 0}

or a dictionary can be used to define them in a model directly. The update_metabolites argument will subsequently update the initial condition in the metabolite object as well.

model.update_initial_conditions({x3: 0, x4:0})
model.initial_conditions

{<MassMetabolite x1 at 0x7fd978e75190>: 1,
 <MassMetabolite x2 at 0x7fd978e75070>: 0,
 <MassMetabolite x3 at 0x7fd978e75340>: 0,
 <MassMetabolite x4 at 0x7fd978e756d0>: 0}

Check the documentation for further details on the MassModel class.

3.4.2. Simulating Dynamic Responses

3.4.2.1. Simulate

Simulating the model once it is set up properly is very simple. To set up the simulation, we use a Simulation object. The simulation object requires a MassModel for initialization.

from mass import Simulation

sim = Simulation(model, verbose=True)

WARNING: No compartments found in model. Therefore creating compartment 'compartment' for entire model.

Successfully loaded MassModel 'Model' into RoadRunner.

The Simulation.simulate method from the will integrate the ordinary differential equations of the system in the provided time interval and return the dynamic responses of concentrations and fluxes.

t0 = 0
tf = 1e6

conc_sol, flux_sol = sim.simulate(
    model, time=(t0, tf), interpolate=True, verbose=True)

Getting time points
Setting output selections
Setting simulation values for 'Model'
Simulating 'Model'
Simulation for 'Model' successful
Adding 'Model' simulation solutions to output
Updating stored solutions

Note: If a model is unable to be simulated, a warning will be raised. By setting the verbose argument to True, a QC/QA report outlining inconsistencies, missing values, and other issues will also be generated and displayed to assist in diagnosing the reason why a model could not be simulated.

3.4.2.2. Inspect the solution

As the default setting, the Simulation object utilizes scipy interpolating functions to capture the concentration and flux responses (see documentation for scipy.interpolate for additional information). The Simulation.simulate_model method returns two cobra.DictLists containing specialized dictionaries known as MassSolution objects.

The first MassSolution object contains the MassMetabolite identifiers as keys, and their corresponding concentration solutions as values.

for metabolite, solution in conc_sol.items():
    print(metabolite, solution)

x1 <scipy.interpolate._interpolate.interp1d object at 0x7fd9495758b0>
x2 <scipy.interpolate._interpolate.interp1d object at 0x7fd9190565e0>
x3 <scipy.interpolate._interpolate.interp1d object at 0x7fd928d44a90>
x4 <scipy.interpolate._interpolate.interp1d object at 0x7fd949583130>

Similarly, the second MassSolution object contains the MassReaction identifiers as keys, and their corresponding flux solutions as values.

for reaction, solution in flux_sol.items():
    print(reaction, solution)

v1 <scipy.interpolate._interpolate.interp1d object at 0x7fd9894511d0>
v2 <scipy.interpolate._interpolate.interp1d object at 0x7fd928d12a90>
v3 <scipy.interpolate._interpolate.interp1d object at 0x7fd989451180>

3.4.2.3. Query time responses

The interpolating functions are functions of time. Therefore, we can evaluate the interpolating function at a specific time point using the following:

time_points = 100;
for metabolite, interpolating_function in conc_sol.items():
    print("%s: %s" % (metabolite, interpolating_function(time_points)))
print()
for reaction, interpolating_function in flux_sol.items():
    print("%s: %s" % (reaction, interpolating_function(time_points)))

x1: 0.3710242389082219
x2: 0.3704524448547136
x3: 0.2569363810507253
x4: 0.0015869284328310024

v1: 0.0005717940535082393
v2: 0.0011351606380398832
v3: 2.5693638105072534e-05

It is also possible to get values for multiple time points at once:

time_points = [0.01, 0.1, 1, 10, 100, 1000];
for metabolite, interpolating_function in conc_sol.items():
    print("%s: %s" % (metabolite, interpolating_function(time_points)))
print()
for reaction, interpolating_function in flux_sol.items():
    print("%s: %s" % (reaction, interpolating_function(time_points)))

x1: [0.99009934 0.90936384 0.56699581 0.4790072  0.37102424 0.32389592]
x2: [0.00990017 0.09058937 0.43018534 0.47682748 0.37045244 0.32388517]
x3: [4.96651050e-07 4.67952432e-05 2.81873928e-03 4.41437817e-02
 2.56936381e-01 3.21735095e-01]
x4: [1.65778860e-13 1.58583691e-10 1.07541682e-07 2.15291648e-05
 1.58692843e-03 3.04838051e-02]

v1: [9.80199167e-01 8.18774471e-01 1.36810476e-01 2.17971609e-03
 5.71794054e-04 1.07505761e-05]
v2: [9.89967148e-05 9.05425714e-04 4.27366598e-03 4.32683701e-03
 1.13516064e-03 2.15007648e-05]
v3: [4.96651050e-11 4.67952432e-09 2.81873928e-07 4.41437817e-06
 2.56936381e-05 3.21735095e-05]

For example, a pandas.Dataframe of concentration values at different time points could be generated using this method:

import pandas as pd

data = [interpolating_function(time_points)
        for interpolating_function in conc_sol.values()]
index_col = [metabolite for metabolite in conc_sol.keys()]
pd.DataFrame(data, index=index_col, columns=time_points)

	0.01	0.10	1.00	10.00	100.00	1000.00
x1	9.900993e-01	9.093638e-01	5.669958e-01	0.479007	0.371024	0.323896
x2	9.900168e-03	9.058937e-02	4.301853e-01	0.476827	0.370452	0.323885
x3	4.966510e-07	4.679524e-05	2.818739e-03	0.044144	0.256936	0.321735
x4	1.657789e-13	1.585837e-10	1.075417e-07	0.000022	0.001587	0.030484

The same can be done for the fluxes:

data = [interpolating_function(time_points)
        for interpolating_function in flux_sol.values()]
index_col = [reaction for reaction in flux_sol.keys()]
pd.DataFrame(data, index=index_col, columns=time_points)

	0.01	0.10	1.00	10.00	100.00	1000.00
v1	9.801992e-01	8.187745e-01	1.368105e-01	0.002180	0.000572	0.000011
v2	9.899671e-05	9.054257e-04	4.273666e-03	0.004327	0.001135	0.000022
v3	4.966510e-11	4.679524e-09	2.818739e-07	0.000004	0.000026	0.000032

3.4.2.4. Filtering for specific species and fluxes

Because concentration and flux MassSolution objects are specialized dictionaries, they can be handled like any other dictionary. Therefore, obtaining the solution for individual species and fluxes can be done easily by using the MassMetabolite or MassReaction identifiers as keys.

print(x1.id, conc_sol[x1.id])

x1 <scipy.interpolate._interpolate.interp1d object at 0x7fd9495758b0>

for flux in [v1, v2]:
    print(flux.id, flux_sol[flux.id])

v1 <scipy.interpolate._interpolate.interp1d object at 0x7fd9894511d0>
v2 <scipy.interpolate._interpolate.interp1d object at 0x7fd928d12a90>

3.4.2.5. Switching between numerical arrays and interpolating functions

Suppose that instead of working with interpolating functions, we would rather work with the original time points and the corresponding solutions utilized by the ODE solver. One way this can be done would be to access the original time point values stored in the Solution object, and use those in the interpolating function:

time_points = conc_sol.t
# Get a slice of the first 50 points
print(conc_sol["x1"](time_points)[:50])

[1.         1.         0.99999943 0.99999734 0.99999525 0.99999316
 0.99998821 0.99997787 0.99995537 0.99990413 0.99978077 0.99942553
 0.99907055 0.99871581 0.99806461 0.99741426 0.99676475 0.9961161
 0.9948862  0.99290131 0.98987466 0.9868666  0.983877   0.98090577
 0.97795277 0.97334374 0.96877916 0.96425858 0.95727034 0.94370569
 0.92186543 0.90109967 0.88135539 0.86258222 0.84473223 0.82775987
 0.81162185 0.78756746 0.76536568 0.7448733  0.72595816 0.70849829
 0.69238106 0.67750262 0.66376717 0.64105858 0.62146897 0.6045664
 0.58997873 0.57738534]

To quickly convert an entire MassSolution object from interpolating functions to numerical arrays or vice-versa, we use the MassSolution.interpolate setter method:

conc_sol.interpolate = False
# Get a slice of the first 50 points
conc_sol["x1"][:50]

array([1.        , 1.        , 0.99999943, 0.99999734, 0.99999525,
       0.99999316, 0.99998821, 0.99997787, 0.99995537, 0.99990413,
       0.99978077, 0.99942553, 0.99907055, 0.99871581, 0.99806461,
       0.99741426, 0.99676475, 0.9961161 , 0.9948862 , 0.99290131,
       0.98987466, 0.9868666 , 0.983877  , 0.98090577, 0.97795277,
       0.97334374, 0.96877916, 0.96425858, 0.95727034, 0.94370569,
       0.92186543, 0.90109967, 0.88135539, 0.86258222, 0.84473223,
       0.82775987, 0.81162185, 0.78756746, 0.76536568, 0.7448733 ,
       0.72595816, 0.70849829, 0.69238106, 0.67750262, 0.66376717,
       0.64105858, 0.62146897, 0.6045664 , 0.58997873, 0.57738534])

conc_sol.interpolate = True
conc_sol["x1"]

<scipy.interpolate._interpolate.interp1d at 0x7fd98944dc20>

for key, value in conc_sol.items():
    print(key, value)

x1 <scipy.interpolate._interpolate.interp1d object at 0x7fd98944dc20>
x2 <scipy.interpolate._interpolate.interp1d object at 0x7fd978e52e50>
x3 <scipy.interpolate._interpolate.interp1d object at 0x7fd978e721d0>
x4 <scipy.interpolate._interpolate.interp1d object at 0x7fd978e72090>

conc_sol["x1"]

<scipy.interpolate._interpolate.interp1d at 0x7fd98944dc20>

conc_sol.x1

<scipy.interpolate._interpolate.interp1d at 0x7fd98944dc20>

3.4.3. Visualizing the Solution Graphically

Once the model has been simulated, the solutions can be visualized using the visualization tools in MASSpy.

import matplotlib.pyplot as plt
import numpy as np

from mass.visualization import (
    plot_phase_portrait, plot_time_profile, plot_tiled_phase_portraits)

All visualization tools utilize the matplotlib python package. See documentation for the visualization class for more details on the available plotting kwargs.

3.4.3.1. Draw time course

Plotting the dynamic responses is straightforward using the plot_time_profile function:

plot_time_profile(conc_sol);

For this model and simulation, plotting on a linear scale does not provide us information about the dyanmics at various time scales. Therefore, we can use the plot_function kwarg to change the scale. Let us keep a linear scale on the y-axis, but change the x-axis to a logarithmic scale.

plot_time_profile(conc_sol, plot_function="semilogx");

The observable argument allows one to specify particular solutions from the solution profile to observe while filtering out all other solutions. For example, only the solutions for \(x_1\) and \(x_2\) can be observed by setting observable to an list of these two keys in the solution profile.

plot_time_profile(conc_sol, observable=["x1", "x2"],
                  plot_function="semilogx");

Though the dynamic behavior is clear, the above plots do not provide any other information. Let us add axes labels, a title, and a legend to the plot.

plot_time_profile(
    conc_sol, legend="right outside", plot_function="semilogx",
    xlabel="Time", ylabel="Concentration",
    title=("Concentration Solutions", {"size": "large"}));

3.4.3.2. Draw phase portraits

Plotting the dynamic responses against one another is also straightforward by using the plot_phase_portrait function:

plot_phase_portrait(conc_sol, x="x1", y="x2",
                    xlabel="x1", ylabel="x2");

\(x_1\) vs \(x_2\): note that you can use the annotate_time_points argument to highlight particular time points of interest. This argument can be utilized either by providing iterable of time points of interest. The annotate_time_points_color can be used to set the color of the time points. To use color to distinguish time points, the number of colors should equal the number of time points specified.

plot_phase_portrait(
    conc_sol, x="x1", y="x2", xlabel="x1", ylabel="x2",
    annotate_time_points=[t0, 1e-1, 1e0, 1e1, 1e3, tf],
    annotate_time_points_color= [
        "red", "green", "purple", "yellow", "cyan", "blue"],
    annotate_time_points_legend="lower outside");

All pairwise phase portraits can be generated and viewed at once in a tiled format using the plot_tiled_phase_portrait function:

plot_tiled_phase_portraits(conc_sol,
                           annotate_time_points_legend="right outside");

This method is particularly useful for looking at correlations at various time scales. For example, looking at the overall behavior, a fast time timescale of (0, 1), an intermediate timescale of (3, 100), and a slow timescale of (300, 10000), we can generate the following:

correlations = [
    np.empty((3, 3)).astype(str),
    np.empty((3, 3)).astype(str),
    np.empty((3, 3)).astype(str),
    np.empty((3, 3)).astype(str)]
for entry in correlations:
    entry.fill("")

fmt_str = "{0:.2f}\n{1:.2f}"
correlations[1][0, 1] = fmt_str.format(*[1, 1])
correlations[2][0, 1] = fmt_str.format(*[1, 1.02])
correlations[2][0, 2] = fmt_str.format(*[1, 0.51])
correlations[2][1, 2] = fmt_str.format(*[1, 0.5])
correlations[3][0, 1] = fmt_str.format(*[1, 1.02])
correlations[3][0, 2] = fmt_str.format(*[1, 1.02])
correlations[3][1, 2] = fmt_str.format(*[1, 1.02])

fig, axes = plt.subplots(2, 2, figsize=(10, 10))
axes = axes.flatten()

times = [(0, 400000), (0, 1), (3, 100), (300, 10000)]
titles = ["{0}\nt0={1}; tf={2}".format(label, *time)
         for label, time in zip(["(a)", "(b)", "(c)", "(d)"], times)]

for i, ax in enumerate(axes.flatten()):
    plot_tiled_phase_portraits(
        conc_sol, observable=["x1", "x2", "x3"], ax=ax,
        plot_tile_placement="lower", additional_data=correlations[i],
        time_vector=np.linspace(*times[i], int(1e6)),
        tile_xlabel_fontdict={"size": "large"},
        tile_ylabel_fontdict={"size": "large"},
        title=titles[i])

3.4.4. Post process the solution

3.4.4.1. Analyze pool behavior

In order to analyze the behavior of pools, pools can be created using the MassSolution.make_aggregate_solution method using the string representation of the pooling formulas. Additional parameters can also be incorporated into the pool formulation using a dictionary input for the parameters argument.

pools = ["x1 - x2", "x1 + x2 - 2*x3", "x1 + x2 + x3"]

for i, equation_str in enumerate(pools):
    pool_id = "p" + str(i + 1)
    conc_sol.make_aggregate_solution(
        pool_id, equation=equation_str, update=True)
    print(pool_id, conc_sol[pool_id])

p1 <scipy.interpolate._interpolate.interp1d object at 0x7fd928d444f0>
p2 <scipy.interpolate._interpolate.interp1d object at 0x7fd9392b5400>
p3 <scipy.interpolate._interpolate.interp1d object at 0x7fd959381310>

This method utilizes the solutions for the individual metabolites over the time range input, and then creates new solutions to represent the behavior of those pools.

plot_time_profile(
    conc_sol, observable=["p1", "p2", "p3"], legend="best",
    plot_function="semilogx",
    xlabel="time",  ylabel="Concentrations",
    title=("Pool profile", {"size": "large"}));

3.4.4.2. Compute and plot the fluxes

A similar process as above can be utilized to obtain behavior of the net flux through a group of reactions. Note that the MassSolution.make_aggregate_solution method relies on the sympy.sympify function and can therefore utilize specific methods, such as the absolute value function, in the string as well.

flux_sol.make_aggregate_solution(
    "v_net", equation='Abs(v1) + Abs(v2) + Abs(v3)', update=True)

{'v_net': <scipy.interpolate._interpolate.interp1d at 0x7fd928ec4f40>}

Again, this method obtains the solutions for the individual fluxes over the time range input, and then creates new solutions to represent the behavior of various flux combinations.

plot_time_profile(
    flux_sol, observable=["v_net"], legend="best",
    plot_function="semilogx",  xlabel="time", ylabel="Fluxes",
    title=("Net Flux", {"size": "large"}));

3.4.4.3. Plot phase portraits of pools

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
plot_tiled_phase_portraits(
    conc_sol, observable=["p1", "p2", "p3"], ax=ax,
    plot_tile_placement="lower",
    annotate_time_points_legend="right outside");

Here, we can see that all of the defined pools are dynamically independent of one another.

3.5. Summary

• Network dynamics are described by dynamic mass balances \((d\textbf{x}/dt = \textbf{Sv}(\textbf{x}; \textbf{k}))\) that are formulated after applying a series of simplifying assumptions

• To simulate the dynamic mass balances we have to specify the numerical values of the kinetic constants \((\textbf{k})\), the initial conditions \((\textbf{x}_0)\), and any fixed boundary fluxes.

• The equations with the initial conditions can be integrated numerically.

• The solution contains numerical values for the concentration variables at discrete time points. The solution is graphically displayed as concentrations over time, or in a phase portrait.

• The solution can be post-processed following its initial analysis to bring out special dynamic features of the network. Such features will be described in more detail in the subsequent notebooks.

\(\tiny{\text{© B. Ø. Palsson 2011;}\ \text{This publication is in copyright.}\\ \text{Subject to statutory exception and to the provisions of relevant collective licensing agreements,}\\ \text{no reproduction of any part may take place without the written permission of Cambridge University Press.}}\)