Modelling

What do we do in computational biological modeling? Simply put, we translate biology into mathematics and back again. It may be to predict, to confirm or to study. Whatever the reason may be, as modelers, our goal is to turn the biology into mathematical terms, which are more easily manipulated, tested and analyzed, and deduce the biological meaning from the results. Above all, our main goal is to work with and support the members conducting the biological research in the laboratory.

The system we designed and built this year is known as a tristable switch (a more detailed description of the biology can be found in the project section). In essence, the system is a two-gene construct with mutual repression and self-activation that should result in three stable states. A diagram of the system is seen in Figure 1. The two genes code for transcriptional factors, GEV and rtTA, which act as self-activators and repressors in the presence of activator molecules (beta-estradiol and aTc, respectively). These small molecules, in effect, act as inducible switches which allow them to control which transcriptional factor is active. Our objectives for the model are to design mathematical models representative of the system in order to 1) predict the stable points of the systems and at what concentrations of beta-estradiol and aTc they occur and 2) to study the dynamics of inducing the switch between stable states. But what does it mean to have stable states?

Simplified visual diagram of genetic construct. Mutual inhibition and self activation are mediated by beta-estradiol and aTc when GEV and rtTa are involved, respectively.

In a biological context, stable states refer to states of stable gene expression levels. This means that expression levels tend to progress until it reaches a certain stable expression level, depending on the initial conditions. These states also are resistant to modifications of expression by changing regulatory proteins, loss of genetic product, etc. In a modeling context, it refers to the trajectory of the system (or group) of equations to converge to a steady or stable state. Mathemtically speaking, stable states occur at intersections of nullclines, but simply, they attract the trajectory of the system of equations. To analyze model stability, one can use phase plane and bifurcation analyses. Phase plane analyses are used to analyze the dynamics of the stability of the system of equations, at a given set of parameters, and bifurcation analyses are used to analyze the change in stability with regards to modifications of parameters.

Phase plane analysis of GEV and rTTA.

In a phase plane analysis, we investigate the stability of two components (molecules like GEV and rtTA) with respect to each other by plotting the amount of the first component (GEV) against the amount of the second (rtTa). For an example, take Figure 2. Each line represents a single simulation, with a unique set of initial conditions (starting amounts) where the blue progressing to red in the line represents the progress of the simulation from the beginning to the end. The phase plane, in essence, is a summary of multiple simulations of the system of equations with different initial conditions. This system of equations tend to progress to three stable states (the groupings of red ends). Stable states are also known as attractor sites or "sinks", because, like a sink, these sites pull in the trajectories of the components as seen by the phase plane.

Bifurcation analysis of GEV and rTTA with respect to parameter a1.

In bifurcation analyses, we investigate the change in stable states with regards to the parameters in order to test robustness (the range of each parameter where the model can still predict tristability) and sensitivity (which parameter can cause the most change in the stable states). The bifurcation diagram is very similar to phase planes; in essence, a phase plane is conducted at each parameter value. At each value, multiple simulations of the system equations are conducted but only the last several values of the simulation (red regions in phase plane) are plotted. One can observe the progress of the stable points as one modifies the parameter. Figure 3 best demonstrates the concept of bifurcation analysis. In this example, we can say that the stable state with high expression of GEV is greatly affected by the parameter.

We just went through why we are modeling and how we're going to use it. The question now is what we are going to use to model? For modeling biological systems, one of the first choices is whether to use a deterministic or stochastic model. In essence, a deterministic model assumes that all variables of the model can account for the majority of the biological behavior of the system, ignoring the inherent variability of the system. Stochastic models take this variability into account by introducing a random factor, but this often makes the system difficult to analyze. We decided to pursue the design of an ordinary differential equation (ODE) deterministic model and a Gillespie-based stochastic model.

Equations and parameters

Deterministic model

Our main objectives with the model are to demonstrate tristability and to study the inducibility of the switch, so we designed ordinary differential equation (ODE) models to achieve the two objectives. ODE models refer to the ordinary differential equation function used in the system which are functions where the derivative of a component is defined with respect to itself (e.g. dx/dt = ax +b, where a and b are constants). In our model, we are restricting it to two equations, each one representing the change in amount of GEV or rtTa with respect to time (dGEV/dt and drtTa/dt), in order to simplify the analysis of stability. In essence, both equations represent the change in expression of the two genes over time and these equations are then solved at each time point to give the expression level at said time point.

Assumptions:

• Ignoring spatial dimensions (nuclear import of protein, nuclear export of mRNA, etc.)
• Grouping of transcription and translation rates into general expression rate
• Ignoring degradation of mRNA

Graphical representation of how each equation was constructed.

For each design of the construct, there were two groups of equations: one for showing tristability and one for studying inducible switch. Within each group, there are three sets of equations, each set handles the terms of self-activation and inhibition differently. Every equation in the ODE models are based on what happens to the promoters driving the expression of GEV and rtTa. There will always be a basal expression rate and maximal expression rate affected by the repression and activation caused by the transcriptional factor.

Schematic diagram of design one of genetic construct

The next three sets of equations represent design one solely for the purpose of demonstrating tristability. The reason for three different sets of equations is to determine which set of assumptions best represents the genetic construct:

Equations for demonstrating tristability of design one assuming multiplicative terms.

Equations (1) and (2) assume that the terms for repression and self-activation are multiplicative. The parameters pertaining to the promoters are: a, which represents the basal expression rate, and b, which represents the maximal expression rate, for pGEVtx in equation (1) and pTREgx in equation (2). The parameters pertaining to the transcriptional factors are k, which represents the degradation rate of proteins, K, which represents the dissociation constant of the transcriptional factor to DNA, and n, which represents the Hill coefficient of binding for the transcriptional factor to DNA. K and n are also specific to the promoter, so the values of n and K for GEV or rtTA are different between equations (1) and (2).

Equations for demonstrating tristability of design one assuming additive terms.

Equations (3) and (4) are similar to equations (1) and (2), except it’s now assumed that the terms of repression and self-activation are additive rather than multiplicative. All parameters are the same to equations (1) and (2), save for d, which represents an arbitrary self-activation strength factor.

Equations for demonstrating tristability of design one assuming singular terms.

Equations (5) and (6) assume that repression and self-activation should be treated as a singular term rather than separate terms. All parameters are the same to equations (1) and (2).

To include the inducible elements of beta-estradiol and aTc to the equations, we assumed that the law of mass action applied and replaced the variables of GEV or rTTA with c_GEV * (GEV * Be) or c_rTTA * (rTTA * aTc), respectively. Be and aTc represent the amount of beta-estradiol and aTc added to the system [molecules] and the parameter c represents the binding constant of the factor and small molecule binding reaction (c = K - 1). Other than that addition, equations (7)-(12) are similar to equations (1)-(6), where equations (7) and (8) assume multiplicative terms, equations (9) and (10) assume additive terms and equations (11) and (12) assume singular terms.

Equations for inducible switch of design one assuming multiplicative terms for design two.

Equations for inducible switch of design one assuming additive terms for design two.

Equations for inducible switch of design one assuming singular terms for design two.

Schematic of design two for genetic construct. Similar to design one, except for the addition of pMRP7 that drives rTTA and GEV expression.

For the second design of the system, there is an addition of the promoter pMRP7 that drives the expression of GEV and rTTA. The addition of a second promoter results in the addition of a second basal expression rate (a2) in each equation from design one. The basal and basal expression rates for pGEVtx (equations (13),(15),(17),(19),(21) and (23)) and for pTREgx (equations (14),(16),(18),(20),(22) and (24)) are now represented by a1 and b1, respectively. This also includes the arbitrary self-activation factor (d1) in equations (15), (16), (21) and (22).

Equations for demonstrating tristability of design two assuming multiplicative terms.

Equations for demonstrating tristability of design two assuming additive terms.

Equations for demonstrating tristability of design two assuming singular terms.

Equations for inducible switch of design two assuming multiplicative terms.

Equations for inducible switch of design two assuming additive terms.

Equations for inducible switch of design two assuming singular terms.

Results