Modeling of the CARTELTM AND Gate

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).

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 (Fig. 1). 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. Therefore 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). Below, it is shown how these ODEs can be set up in their simplest form.

$f\text{'}\left(t\right)=a\left(z\right)-{k}_{1}*f\left(t\right)$

(1)

$y\text{'}\left(t\right)={k}_{2}*f\left(t\right)-{k}_{3}*y\left(t\right)$

(2)

The variables are functions of time $t$, where $f\left(t\right)$ describes the mRNA concentration, $a\left(z\right)$ the transcription function, ${k}_{1}$ the mRNA degradation rate, $\mathrm{y\left(t\right)}$ the protein concentration, ${k}_{2}$ translation rate and ${k}_{3}$ 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 the delay of protein expression is not relevant for this model.

$x\text{'}\left(t\right)={k}_{4}+{k}_{5}*b\left(t\right)-{k}_{6}*x\left(t\right)$

(3)

The variables are functions of time $t$, where $x\left(t\right)$ describes the protein concentration over time, ${k}_{4}$ the basal expression rate, ${k}_{5}$ the maximum expression rate, $b\left(t\right)$ the activation function and ${k}_{6}$ the degradation rate.

Regarding our intention to raise the specificity of cell killing, different possible designs were conceivable. These are accompanied by the introduction of certain genes or their elimination through a knockout. Focusing on the controlled CAR expression, an AND gate with a clear defined allocation of inputs and output had to be generated.

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 promoter. 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 VEGFR-2. 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 (equations in following drop down menus, Fig. 2, 3).

${k}_{basal}$ describes the basal expression rate, ${k}_{max}$ the maximum expression rate, ${k}_{deg}$ 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. ${k}_{dim}$ and ${k}_{dis}$ are always describing the dimerization and dissociation rate of HIF1AB.

${k}_{basal}$ describes the basal expression rate, ${k}_{max}$ the maximum expression rate, ${k}_{deg}$ 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. ${k}_{dim}$ and ${k}_{dis}$ are always describing the dimerization and dissociation rate of HIF1AB.

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 (Fig. 4). Therefore the HIF1A knockout was chosen.

After deciding for an AND gate design, we got in contact with the group of Professor Jens Timmer. Cooperating with Svenja Kemmer from his group, we gained insights into the methodology of processing and evaluating experimental data by modeling. In order to model our system of interest, the previously found 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. The software used for the modeling of ODEs is data2dynamics, which is applied in Matlab and has been developed by the group of Professor Jens Timmer.

The composition of the model is described below in Figure 5.

The respective equations are equivalent to those of the first model except for an additional equation describing the HIF1A mRNA since the CRE characterization data showed delay in the processing.

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 system.

Additionally, the dependency of the HIF1A concentration on the pH and the oxygen concentration is shown in a heat map (Fig. 6). 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. Since 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.

The mathematical description of our AND gate leads to more ideas of how to improve the design of our logic gate. For example an improvement of the AND gate could be achieved by varying the number of enhancer elements. Depending on whether there is any cooperation in the binding to the enhancer elements and 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.

Effective Range of Activated CARTELTM T cells

We often received the question 'How far can a T cell migrate after being activated in a tumor microenvironment?'. We consulted this year’s iGEM team from ETH Zurich due to their modeling expertise to find an answer to this issue. Their basic model indicated that cells migrate on the average 6 mm from the tumour microenvironment before getting inactivated.