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 parameters to generate a desired output if enough data is available (Chen et al., 1999).

Figure 1: Schematic description

Finding the CARTELTM AND gate

The hypoxia response element (HRE), cAMP response element (CRE) and CTLA4 promoter are responsive to traits specific for the tumor microenvironment. They can therefore be utilized for a more specific CAR expression. Combining these three enhancers using an AND gate to allow a specific CAR expression (Brophy et al., 2014) proved to be quite challenging. As T-cells natively use and modulate these sensing mechanisms crosstalk would have been inevitable. A knockout or knockdown of either HIF1, VEGFR-2 or TDAG8 had to be generated to avoid interference from endogenous signaling.

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 randomly chosen AND gate design.

To create such a model, the changing concentration of CAR can be described by a system of coupled ODEs. For each component of the genetic circuit - proteins and mRNAs - rate equations describing their respective 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 provided or 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 as it can be described as 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 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]

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 knockout possibilities In the following (1) stands for high and (0) for low induction level. The first number represents the induction level via VEGF or pH and the second via oxygen concentration. The colors symbolize the following conditions: Red (1/1), Blue (1/0), Yellow (0/1), Green (0/0). a) Depiction of relative CAR expression under all four conditions over the curse of 24 h in a HIFA knockout. It was assumed that at 0 h relative CAR expression is at 0 %. b) Depiction of relative CAR expression under all four conditions over the curse of 24 h in a HIFA knockout. It was assumed that at 0 h relative CAR expression is at 100 %. c) Depiction of relative CAR expression under all four conditions over the curse of 24 h in a TDAG8 of VEGFR-2 knockout. It was assumed that at 0 h relative CAR expression is at 0 %. d) Depiction of relative CAR expression under all four conditions over the curse of 24 h in a TDAG8 of VEGFR-2 knockout. It was assumed that at 0 h relative CAR expression is at 100 %.

The in silico induction of the different knockouts possibilities revealed that a TDAG8 or VEGFR-2 knockout would have a higher relative CAR expression under the conditions of high VEGF and low oxygen concentration than a HIF1A knockout. Therefore the HIF1A knockout was chosen.

After deciding on 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 a 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 relative amounts 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 behavior.

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.