Team:LCG-UNAM-Mexico/Simulation

 
 

 
With the aim of predicting the behavior of the system, the biochemical reactions were implemented in the SimBiology package of MATLAB, using the previously defined parameters. Simulations were run for different values of the initial concentration of AHL and Ni_total (Ni_int + Ni_ext) which are the metabolites that we can directly manipulate in our experiments. A parameter scan was also run for some parameters to understand their influence on the system.
In order to gain insights into the system dynamics to elucidate the conditions needed to get the desired behavior, we performed a series of analyses on it: Sensitivity analysis allowed us to identify critical parameters that needed to be defined on the most stringent way. Basis for the (right) null and left null space were calculated to obtain information about the general network behavior. Steady-states were calculated by numerical integration of the non-linear ODEs system. Finally the Jacobian of the system was calculated around the steady-states. All simulations and analyses were implemented and performed on MATLAB.

Simulation and parameter scan

Once the model was implements on SimBiology, simulations were carried out to analyze its behavior through time after a given AHL concentration. We would like to observe a moderate and ephemeral change in the intracellular nickel concentration; this would allow our bacteria to sing without killing them.

After exhausting our literature search for parameters, we tuned the remaining ones in such a way that they ensure the desired behavior. We performed a parameter scan on the CI dimerization constant (an early step in the signaling cascade) and the rates at which nickel enters and leaves the cells.

CI dimerization

Several scans where performed on the CI rates of dimerization and dissociation to define the range in which the desired response would be observed.

[CI] equilibrium is defined by its synthesis and degradation rate, while k_4.1/k_-4.1 determine the equilibrium of [CI:CI] and modify time that takes [CI] to reach a steady-state.

Several wide range scans show that the behavior of metabolites concentrations is determined by the ratio k_4.1/k_-4.1 and not their specific values; the bigger this ratio the bigger the final dimer concentration, which turns into a greater repression of rcnA and nickel internalization; the time that takes the system to reach a new equilibrium after being perturbed by AHL also increases with this ratio.

We tested the effect of AHL on CI and its dimer. The graph shows that [AHL] only defines the strength of the response without affecting the equilibrium of the CI molecules.

A constant repression of rcnA would be toxic for bacteria, for this reason k_4.1/k_-4.1 cannot be too big. We choose the values k_4.1=0.0001 molecules^-1s^-1 and k_-4.1=0.01s^-1 for their similarity with other dimerization parameters (reaction 2.1) and the fact that they show the desired behavior. Although it must be made clear that there is a wide range of values that would work for our purposes.

These values are comparable to others typical biochemical parameters. It has been shown that kinetic parameters can be modified by changing amino acid sequences (for example, CI half life is reduced by adding a LVA tail in the C-terminal), it’s proposed that it’s possible to engineer the protein to reach an acceptable dissociation constant.

Nickel internalization and efflux

We performed several scans on the rates of nickel internalization and efflux to define the range of values that show the desired response.

Just like in the previous case, our results indicate that the behavior is determined by the the ratio k₇/k₉ and not to their specific values. The smaller the ratio greater the internal nickel concentration.

[Unk] is also an unknown parameter, but as we assume it constant through time it can be completely arbitrary, as long as the reaction rate k₉ is well defined. For convenience we define [Unk] as one tenth of initial [RcnA] due the fact that in the real system we will have 10 more copies of rcnA than in wildtype cells.

We define k₇=500 molecules^-1s^-1 which is a known* extrusion rate, and k₉=500 molecules^-1s^-1 which show the desired behavior. Nonetheless the range of parameters that fit our interests is very large and we think that the real parameters will be within it.

Simulation

Once all the parameters were defined, we simulated the model for 20,000 seconds using a variable order solver based on the numerical differentiation formulas, ode15s(Referencia 2 de mariana).

We performed a scan on the initial AHL concentration to analyze the response of the system.

It can be seen that the final response depends on our input signal and we can manipulate its strength and length by modifying the initial AHL concentration.

A scan over the total nickel concentration shows that response is only weakly influenced by it; even though the intensity of the response does vary, the proportion with respect to total nickel does not change. We decided to fix the initial value on 30,000 molecules for the following analysis, but this number will be finally determined by the bacterial tolerance and the sensitivity of our measurement device.

Sensitivity analysis

A sensitivity analysis using Simbiology was performed to determine the parameters that influenced our model the most. This analysis was possible because our model was solved using an ODE solver.  The values for the calculated sensitivities were fully normalized to dedimensionalize them. All of the parameters were plotted vs. all of the species and the sensitivity at each intersection was determined. There were two main evaluations, one at t=7700s which is approximately the time at which the greatest response was observed.

Sensitivity Analysis, at a time of 7700, at which our model predicts the greatest response.

At that time, the species influenced the most are AHL, cI, cI:cI, RcnA and Ni-ext. That makes sense for it is at this point of time where the production of cI, its dimerization followed by the repression of RcnA and therefore the decrease in Ni-ext, are at the highest levels.

The other analysis at t=74000s which is approximately the time the system takes to return to a basal state.

Sensitivity Analysis, at a time of 74000, at which our model predicts a return to the basal state.

At this time the system has returned to the basal state. Here the species influenced the most are the same; the difference resides in the magnitude of the influence.

Afterwards the sensitivity of a single specie with respect to all of the parameters was plotted with respect to time.

 Time Dependent Sensitivity Analysis, specific for AHL with respect to all of the parameters

The parameters that influenced AHL the most were those involved on its degradation and the formation of the complex AHL:LuxR.

Time Dependent Sensitivity Analysis, specific for cI:cI with respect to all of the parameters

An interesting behavior can be seen for the dimer cI:cI. It is sensitive to most of the parameters in our model.  However most of them are only temporal, only influencing the model during the lifetime of AHL. The only parameters that always influence its concentration are its dimerization and degradation rate and the “leaky” synthesis rate.

Time Dependent Sensitivity Analysis, specific for RcnA with respect to all of the parameters

 RcnA is sensitive to practically the same parameters as cI:cI, nevertheless, it is easily noticeable that RcnA is inversely affected by them, which makes sense, because  the dimerization of cI lowers the available amount of cI and enhances the synthesis of RcnA

 
Time Dependent Sensitivity Analysis, specific for intracellular nickel with respect to all of the parameters

In our model, an increase in the amount of internal nickel is the expected outcome after the addition of AHL. The sensitivity analysis demonstrates that it is sensitive to all of the parameters in our model. Just temporally to some of them, but continuously to the rest.

Stoichiometric matrix

Mathematically, the stoichiometric matrix S is a linear transformation of the flux vector to a vector of the time derivatives of the concentration vector of the form dx/dt = Sv, and it’s formed from the stoichiometric coefficients of the reactions that comprise a reaction network. This matrix is organized such that every column corresponds to a reaction and every row corresponds to a compound. Each column that describes a reaction is constrained by the rules of chemistry, such as elemental balancing; and every row thus describes the reactions in which that compound participates and therefore how the reactions are interconnected. The stoichiometric matrix thus contains chemical and network information.

There are four fundamental subspaces associated with a matrix, and each of them has important roles in the analysis of biochemical reaction networks. In particular, we have analyzed the null space and the left null space of the stoichiometric matrix of the system (Suppl. 1). The null space of S contains all the steady-state flux distributions allowable in the network. The left null space of S contains all the conservation relationships, or time invariants, that a network contains.1(Palsson)

The null space of S was calculated for the system (Suppl. 1). It shows that to reach the steady-state of the system we need to meet the following conditions:

v₄=v_3.1+v_3.2, which means that the total synthesis of CI should be equal to its degradation;
v₆=v₈, the synthesis of RcnA equal to its degradation; and
v₇=v₉, the rate of nickel uptake equal to its efflux rate.

The left null space of the systems shows 4 metabolites with constant concentration, AiiA, ρ, ρCI and Unk, and two mass conservation cycles (moieties), Ni_int+Ni_ext and LuxR + AHL:LuxR + 2((AHL:LuxR):(AHL:LuxR)).

Steady states

The steady state of a system is defined as a nonequilibrium state through which matter is flowing and in which all components remain at a constant concentration^ref. We defined the time taken to reach the steady state as the time when the variation in the concentrations of all components cease to be biologically relevant (here defined as variation in less than one molecule). We can have two initial states in our system:

With no AHL. In this case, the only reaction that takes place is the uptake and efflux of nickel, which instantaneously reaches the steady state.
With AHL. Here, all the regulatory steps in the cascade take place and the steady state is reached when all AHL has been degraded and the signal has faded.

In both cases, the steady state reached is the same. The reason for this is that the system is a regulatory cascade and does not have any feedback loops, so it is expected to have only one steady state. Here we present the simulation in 100,000 seconds, you can click on in to see a larger version:

The concentrations at which it is reached are the following:

The range of AHL concentration in this analysis was chosen from 30 (when we start to see a measurable response) to 10,000. We can observe a linear dependence with higher AHL concentrations. The extent of RcnA repression can be adjusted at will as long as the extracellular concentration of nickel is not damaging to the cell.

Jacobian

The Jacobian of a system is defined as the matrix of first order partial derivatives, ad it represents the best linear approximation to a function at a given point. In biochemical networks, the Jacobian can be defined for metabolites(Jx) or fluxes(Jv):

Where S is the stoichiometric matrix and G is the gradient matrix. S defines the structure of the network and has the stoichiometric coefficients of all reactions which are represented by the rows of S while metabolites are represented by columns. G is the matrix of first order derivatives of fluxes with respect to species concentrations:

This formal representation of the Jacobian formalizes the relation between the topology of the network, and its biophysical and kinetic characteristics (REFERENCIAS DE SUR UNO Y DOS).

For our system, we first obtained S which is rank deficient: its rank is 7 and it has 13 rows, this difference is explained by the moieties of the system. We reduced S to make it congruent with its rank by eliminating the rows corresponding to: AiiA, LuxR, pcI, p, Ni-ext and Unk. Then we constructed G as previously defined; the partial derivatives are calculated only with respect to the species that remained in S; the state at which G is calculated is the steady-state.

Both Jacobians were calculated and they were decomposed through similarity transformation of their eigenvalues and eigenvectos to obtain the modal matrices. This matrices represent the formation of molecular ‘pools’ on the system and the negative inverse of the eigenvalus represent the time scales at which this interactions occur. The Jacobians and modal matrices can be downloaded on Excel (.xls) format if you click here.

The values in the tables show that the nickel response is very fast, in the order of 6.03e-8 s, which is what we need to make this system and efficient transcriptional indicator. On the other hand, there is a pool of RcnA that forms immediately and another one that takes much more time to form.

REFERENCIAS DE SUR
Jamshidi N and Palsson (2008) Formulating genome-scale kinetic models in the post-genome era Molecular Systems Biology 4:171.
Jamshidi N and Palsson (2008) Top-Down analysis of temporal hierarchy in biochemical reactions networks PLoS Comput Biol 4(9).

REFERENCIAS DE MARIANA
1. *Referencia Dan. 2. SimBiology; User’s guide.