Team:Freiburg/Model


Modeling

In synthetic biology, modeling can be applied to a wide range of topics, for example modeling of genetic circuits. In this subfield of mathematical and computational biology, ordinary differential equations (ODEs) are used to describe the transcriptional and translational processes over time, predict the behavior of the desired circuit and also to support their further development such as optimizing for a desired output if enough data is available (Chen et al., 1999).

Figure 1: Schematic depiction

Finding the CARTELTM AND gate

The hypoxia response element (HRE), cAMP response element (CRE) and CTLA4 promoter can be utilized for a more specific CAR expression. Combining these three enhancers in an AND gate to allow a specific CAR expression (Brophy et al., 2014) proved to be quite challenging. Using endogenous systems came with a drawback - a knockout or knockdown of either HIF1, VEGFR-2 or TDAG8 had to be generated.

Generating a single knockout or knockdown is challenging enough but generating three to test all possibilities seemed quite impossible in six months. So using modeling to find the best combination of enhancers and knockout was an attractive option instead of choosing a random AND gate design.

To create such a model, the changing concentration of CAR can be described with a system of coupled ODEs. For each component of the genetic circuit - proteins and mRNAs - rate equations describing its synthesis, degradation, association or other processes are added (Chen et al., 1999). How these ODEs can be set up in their simplest form is shown below.

f'(t)=a(z)-k1*f(t)

(1)

y'(t)=k2*f(t)-k3*y(t)

(2)

The variables are functions of time t, where f(t) describes the mRNA concentration,a(z)the transcription function, k1 the mRNA degradation rate, y(t)the protein concentration, k2 translation rate and k3 protein degradation rate. This kind of ODEs can be solved via numerical integration to describe changing protein and mRNA levels. To obtain the rate constants, experimental data have to be produced which is suitable for a nonlinear regression (Ingalls et al., 2012).

Equations (1) and (2) were further simplified for the AND gate comparison into equation (3) by neglecting the mRNA level because it better describes delayed protein expression.

x'(t)=k4+k5*b(t)-k6*x(t)

(3)

The variables are functions of time t, where x(t) describes the protein concentration over time, k4 the basal expression rate, k5 the maximum expression rate, b(t) the activation function and k6 the degradation rate.

The variables are functions of time t, where x(t) describes the protein concentration over time, k4 the basal expression rate, k5 the maximum expression rate, b(t)the activation function and k6 the degradation rate.

The first possible system was a knockout of HIF1A or HIF1B and its reintroduction via lentiviral transduction under the control of either CRE or CTLA4. Considering the law of mass action, a knockout of HIF1B would have a kinetic disadvantage over a knockout of HIF1A because HIF1B is accumulated permanently and independently from external conditions (Pescador et al., 2005) while HIF1A is only accumulated under hypoxia. Therefore only a knockout of HIF1A was considered a possible candidate for an AND gate.

The second or third possibility would be a knockout of either TDAG8 or VEGFR2. Both can be described by the same model as both receptors are activated by one input and trigger activation of the associated promoter afterwards.

Unfortunately, experimental data only exist for HEK cells which are not representative for T cells (Ausländer et al., 2014; Nguyen et al., 2012). Therefore, the different AND gate designs were modeled with the same rate constants to compare their performances.

Figure 2: Schematic representation of an HIFA knockout AND gate utilizing either low pH or high VEGF concentration as first input and low O2 concentration as second input to express the CAR

kbasal describes the basal expression rate, kmax the maximum expression rate, kdeg the degradation rate, Km gives the concentration of the activator were the promoters half activity is reached and the Hill coefficient n the activational slope of the given promoter. kdim and kdis always describes the dimerization and dissociation rate of HIFAB.

d[HIF1A](t)dt=kbasal+kmax*[VEGF/H+]nKmn+[VEGF/H+]n-kdeg*[HIF1A]

-[HIF1A]* f([O2])- kdim*[HIF1A]*[HIF1B] + kdis*[HIF1AB]

d[HIF1B](t)dt=kdis*[HIF1AB] - kdim*[HIF1A]*[HIF1B]

d[HIF1AB](t)dt =kdim*[HIF1A]*[HIF1B] - kdis*[HIF1AB]

d[CAR](t)dt =  kbasal+ kmax*[HIF1AB]nKmn+[HIF1AB]n - kdeg*[CAR]



Figure 3: Schematic representation of an VEGFR-2 or TDAG8 knockout AND gate utilizing low oxygen concentration as first input and either low pH or high VEGF concentration as a second input to express the CAR

kbasal describes the basal expression rate, kmax the maximum expression rate, kdeg the degradation rate, Km gives the concentration of the activator were the promoters half activity is reached and the Hill coefficient n the activational slope of the given promoter. kdim and kdisalways describes the dimerization and dissociation rate of HIFAB

d[HIFA](t)dt=kbasal-kdeg*[HIF1A]-f([O2])*[HIF1A]

d[HIF1B](t)dt=kdis*[HIF1AB] - kdim*[HIF1A]*[HIF1B]

(10)

d[HIF1AB](t)dt =kdim*[HIF1A]*[HIF1B] - kdis*[HIF1AB]

d[TDAG8/VEGFR2](t)dt=kbasa+kmax*[HIF1AB]nKmn+[HIF1AB]n-kdeg*[TDAG8/VEGFR2]

d[CAR](t)dt=kbasa+kmax*([TDAG8/VEGFR2]*[VEGF/H+])nKmn+([TDAG8/VEGFR2]*[VEGF/H+])n-kdeg*[CAR]


Figure 4: Modeling of knockouts The modeling of different knockouts showed leakiness for receptor knockouts (Fig. 4) whereas the HIF1A knockout has a better selectivity. Therefore, we decided to combine the AND gate with a knockout of HIF1A.

After deciding for an AND gate design, we got in contact with the group of Professor Timmer. Cooperating with Svenja Kemmer from his group, we gained insights into the methodology of processing and evaluating experimental data by ODE modeling. In order to model our system of interest, data of CRE and HRE characterizations in HEK cells were used (Ausländer et al., 2014; Nguyen et al., 2012). Subsequently, parameters as kinetic rates were estimated by the Maximum Likelihood method and used to further analyze the system.

Cell Culture

Figure 5: Representation of the steps describing the input processing in the AND gate

The composition of the model is described in Figure 5.

The respective equations are equivalent to those of the first model except of an additional HIF1A mRNA. As the CRE characterization data showed delay in the processing, this step was added. The mathematical model is able to describe dynamics as well as dose response data. However, experimental measurements were only available for a few components of the system. Thus, detailed predictions require further investigations. Performing promoter characterizations and measuring degradations rates will enable more reliable predictions.

The respective response curves were fitted and combined in the AND gate model, thus describing the functionality of our own system. In the following, the dependency of the HIF1A concentration on the pH and the oxygen concentration is shown in a heat map. All data refer to a time of 24 hours after exposure to the respective conditions. The HIF1A concentration depends on both a low pH and hypoxia as it is present in the tumor microenvironment. As it regulates the CAR expression, it can be estimated that the promoter for the CAR expression will only be activated if both conditions are fulfilled.

Heat map

Figure 6: Heatmap
Illustration of the model predictions of the HIF1A concentration in dependence of the pH and the oxygen concentration in the microenvironment. Different colors describe the concentrations of HIF1A specified in the sidebar.

The mathematical description of our system leads to more ideas of how to improve the design of our logic gate. A possible improvement of the AND gate could be achieved by varying the number of enhancer elements. Given any effect of the cooperation of enhancer elements or the leakiness of the promoter, different constructions are conceivable. At this point, modeling becomes even more important to create an encompassing overview about the different possibilities and simulate the system’s behaviour.

If you are interested in trying out different functions of data2dynamics applied to our AND gate model, you can download the code for it here. To run the model, also Matlab is required.