In order to enable an in-depth understanding of the mechanisms underlying our system, we developed a detailed mathematical model. Our system is based on quorum sensing: Our engineered CAR T-cells constantly express a ligand, but only at a low basal transcription rate. Upon CAR T-cell clustering however, the ligand concentration in the cellular environment exceeds a threshold, so that a positive feedback loop is activated, increasing the ligand and - importantly - the CAR transcription rate. Figure 1 provides a detailed description of our system, with a focus on what is important in order to derive a mathematical model.

Figure 1: Schematic representation of our quorum sensing loop system.

1. MESA dimerization: A ligand triggers two nanobodies to dimerize. As a result, the transcription factor tTA is cleaved by a protease.
2. tTA dimerization: Two cleaved tTA dimerize.
3. DNA binding: The dimerized tTA binds to the Tet Response Element (TRE). This activates ligand and CAR expression.
4a. Ligand: The ligand is expressed and secreted.
4b. CAR: The Chimeric Antigen Receptor is expressed and transported to the cell membrane.

Figure 1: Schematic representation of our quorum sensing loop system.

1. MESA dimerization: A ligand triggers two nanobodies to dimerize. As a result, the transcription factor tTA is cleaved by a protease. 2. tTA dimerization: Two cleaved tTA dimerize. 3. DNA binding: The dimerized tTA binds to the Tet Response Element (TRE). This activates ligand and CAR expression. 4a. Ligand: The ligand is expressed and secreted. 4b. CAR: The Chimeric Antigen Receptor is expressed and transported to the cell membrane.

The following list provides an overview of the proteins important for the functionality of our quorum sensing feedback loop:

• \(M\): Dimerized MESA receptor on the cell surface, consisting of GFP and/or mCherry nanobodies
• \(L\): Ligand, consisting of two domains: GFP-mCherry and/or 2xmCherry
• \(T\): Transcription factor (TF) tTA
• \(C\): Chimeric Antigen Receptor (CAR)

Following the approach of Kannan, 2018 [1] , we describe the general rate equation of a protein with concentration \(x\) as follows: \[\frac{dx}{dt} = synthesis \pm transformation \pm transport − (degradation + dilution)\] The transcription and translation of the protein are included in the synthesis, the transformation accounts for processes like the dimerization of the protein, and the transport considers the molecule diffusion within the cell or through the cell membrane. Moreover, the formula takes into account that the protein concentration might decrease due to degradation, cell dispersion, and dilution of cell culture.

Based on this generic expression, we derived an individual formula for the general rate equation of each protein that is part of our quorum sensing system. For this, we will use the symbols in the following table:

In the following section, we provide a rate equation for the concentration of each protein in our loop system. We assume the MESA receptor \(M\) to be transcribed at a constant rate \(\mu_{bM}\). After expression, the receptors are transported to the cell membrane to be presented there with transportation rate \(\theta_M \). Additionally, we take into account a protein degradation at first-order degradation rate \(\delta_M\). \[\frac{d[M]}{dt} = \mu_{bM} \times \theta_M − \delta_M \times [M] \]

Based on the rate equation of \([M]\), we derive the rate equation for the concentration of dimerized MESA receptor \(M_{di}\) on the cell surface. The dimerization occurs upon binding of an extracellular ligand and is therefore additionally dependent on \([L_e]\). The rate constant for the MESA-ligand complex formation is denoted as \(K_{ML_e}\). \[\frac{d[M_{di}]}{dt} = K_{ML_e} \times [M]^2 \times [L_e] \]

We distinguish between intracellular ligand concentration \([L_i]\) and extracellular ligand concentration \([L_e]\), and assume the ligand is expressed at a basal transcription rate \(\mu_{bL}\). Upon binding of the transcription factor tTA to the DNA binding motif TER, the ligand expression is enhanced. Following the approach of Kannan, 2018 [1], we assume that the TF-DNA binding process follows Michaelis-Menten kinetics with \(K_{T}\) describing the equilibrium constant associated with the binding. We take into account a ligand secretion at rate \(\lambda_L\) and a first-order degradation at rate \(\delta_L\). \[\frac{d[L_i]}{dt} = \mu_{bL_i} +\mu_{mL_i} \frac{[T_{di}]}{k_{T_{di}} + [T_{di}]} - \lambda_L [L_i] - \delta_L [L_i] \]

For the extracellular ligand concentration, we take into account the number of cells \(N\), which secrete ligands at secretion rate \(\lambda_L\). Furthermore, we assume a decrease in extracellular ligand concentration due to ligands that diffuse away from the cell cluster at the diffusion rate \(\phi_{Laway}\). \[\frac{d[L_e]}{dt} = N \times \lambda_L [L_i] - \phi_{Laway}[L_e] \]

Upon dimerization of two nanobodies, resulting in \(M_{di}\), one tTA monomer gets cleaved by the protease. We assume two tTA monomers to dimerize at rate \(K_{diT}\). We account for tTA degradation at first-order degradation rate \(\delta_{T}\). \[\frac{d[T_{di}]}{dt} = [M_{di}]^2 \times K_{diT} - \delta_{T_{di}} [T_{di}]\]

Upon Binding of tTA to the DNA binding site TER, the ligand and CAR transcriptions are enhanced. As already described, we assume that the TF-DNA binding process follows Michaelis-Menten kinetics, with \(k_{T_{di}}\) describing the equilibrium constant. After expression, the CAR is transported to the cell membrane to be presented there with transportation rate \(\theta_C\). CAR degradation happens at first-order degradation rate \(\delta_C\). \[\frac{d[C]}{dt} = (\mu_{bC} +\mu_{mC} \frac{[T_{di}]}{k_{T_{di}} + [T_{di}]}) \times \theta_C - \delta_C [C] \]

In the following, we will use the system of ordinary differential equations derived in the previous section to mathematically simulate our quorum sensing feedback loop. The goal is to obtain a numerical simulation that depicts the CAR concentration on the cell surface \([C]\) with respect to the cell number \(N\).

We solved the equations using the nagf_ode_ivp_bdf_zero_simple (d02ejf) method of the python interface of the Numerical Algorithms Group (NAG), which implements backward differentiation. This procedure is in accordance with the approach of Ward et al., 2001 [2].

We assumed a starting concentration of 0.002 for all our proteins except for the not dimerized MESA receptors \(M\), since the nanobodies are expressed irrespective of the cell number \(N\). We estimated a starting concentration of 0.65 for \(M\). Furthermore, we performed literature research to obtain reasonable values for the rate constants. The parameters we used for our numerical simulation can be found below.

Results

Figure 2: Numerical simulation of our quorum sensing loop system.

The magnitude of the results, especially with focus on the cell number, matches the findings of Kannan, 2018, Figure 5 [1], who also simulated a system based on quorum sensing and implemented the model with similar parameter values.

However, as most of our parameters are estimations, our numerical simulation has to be interpreted as a rough approximation that does not reflect the complexity of all simultaneous processes within a cell cluster.

Our numerical simulation of the CAR concentration [\(C\)] depicts the expected behavior: the effect of our feedback loop sets in only if several cells cluster together. The exact threshold for the feedback loop activation is highly dependent on the estimated parameters. Consequently, the threshold can be altered by fine-tuning those rate constants.

In the following sections, we will perform an in-depth investigation of the docking processes determining the parameters \(k_{T_{di}}\) and \(K_{ML_e}\). Furthermore, we use mutagenesis to alter the docking affinity, which allows for precise fine-tuning of our loop activation threshold.

In the following section, we examine the binding of our dimerized transcription factor (TF) tTA to the DNA binding motif TRE (Figure 1, step 3). The TF tTA is composed of a tetR DNA binding domain and a VP16 transcription activation domain [3]. That way, the transcription of the controlled gene (in our case the genes encoding the ligand and the CAR) is enhanced upon tTA-DNA binding, whereas the binding of the tetR domain without the VP16 activation domain would repress the transcription.

The TF-DNA binding affinity is represented by the rate constant \(k_{T_{di}}\) in our mathematical model. As described before, the loop system activation threshold can be altered by optimizing the equilibrium constant of this binding process, for instance by making the binding more or less energetically preferable. This fine-tuning can be achieved by energetically assessing the binding of tTA and the DNA binding motif TRE and subsequently mutating the binding residues of the TF and the bases of the TRE, while investigating how this changes the binding affinity. Of note, all considerations in this section assume the absence of tetracycline. An investigation of the tetracycline-TF docking can be found in the subsequent section.

Code availability: The code we wrote to perform the mutation of the TF and the molecular docking can be found on our Software GitLab.

Handling of signal sequences: Some of our sequences encoding the proteins relevant for our loop system contain a signal sequence, which, for example, results in the subcellular translocation of the protein. As this is irrelevant to the binding process, the analysis on this website is restricted to protein sequences without signal sequences. All sequences provided on this website exclude the signal sequence.

Obtaining a PDB structure for the transcription factor

As portrayed in Figure 1, step 2, the TF tTA binds as a dimer to the TRE DNA binding motif [4]. The first step towards simulating the TF-DNA binding kinetics is to obtain a 3D structure of the dimerized TF.

As described, the TF tTA is composed of a tetR DNA binding domain and a VP16 activation domain. As the VP16 activation domain is irrelevant to the TF-DNA binding kinetics, and has an unknown structure, we restricted our further downstream analysis to the TetR domain of the TF.

▼ View tTA sequence

MSRLDKSKVINSALELLNEVGIEGLTTRKLAQKLGVEQPTLYWHVKNKRALLDALAIEMLDRHHTHFCPLEGESWQDFLRNNAKSFRCALLSHRDGAKVHLGTRPTEKQYETLENQLAFLCQQGFSLENALYALSAVGHFTLGCVLEDQEHQVAKEERETPTTDSMPPLLRQAIELFDHQGAEPAFLFGLELIICGLEKQLKCESGSAYSRARTKNNYGSTIEGLLDLPDDDAPEEAGLAAPRLSFLPAGHTRRLSTAPPTDVSLGDELHLDGEDVAMAHADALDDFDLDMLGDGDSPGPGFTPHDSAPYGALDMADFEFEQMFTDALGIDEYGG

▼ View TetR part of the tTA sequence

MSRLDKSKVINSALELLNEVGIEGLTTRKLAQKLGVEQPTLYWHVKNKRALLDALAIEMLDRHHTHFCPLEGESWQDFLRNNAKSFRCALLSHRDGAKVHLGTRPTEKQYETLENQLAFLCQQGFSLENALYALSAVGHFTLGCVLEDQEHQVAKEERETPTTDSMPPLLRQAIELFDHQGAEPAFLFGLELIICGLEKQLKC

The 3D structure of the TetR protein has been experimentally confirmed: The Protein Data Bank (PDB) entry 4AC0 has 100% sequence identity with our TetR sequence. We therefore took the PDB file of entry 4AC0 and subsequently obtained a 3D structure of the dimerized TetR using the HDOCK server, which is based on a hybrid algorithm of template-based modeling and ab initio free docking [5].

Figure 3: Result of the TetR dimerization using HDOCK. The video was generated with PyMol. Tool: HDOCK; Input Receptor Molecule: 4AC0.pdb; Symmetric multimer docking: C2.

Identifying the TetR binding residues and the DNA binding site

Using PredictProtein [6] to analyze the sequence of the TetR domain of our TF, we identified the residues 39-44 as primarily responsible for the DNA binding.

▼ View PredictProtein results

The DNA binding motif of TetR is known [4]. We identified the binding motif in our DNA sequence and centered the binding motif within a sequence of 29 bases:

ACCACTC CCTATCAGTGATAGA GAAAAGT

We started with investigating the binding affinity between the dimerized TetR protein and the DNA binding motif, both without any mutations. We performed the docking simulation with HDOCK, using the following settings:

• Input: pdb file of the dimerized TetR protein
• Receptor binding site residues: 36-47:A, 36-47:B
• Ligand binding site residues: 8-11:A, 17-19:A, 37-40:A, 46-48:A

The four ligand (i.e. DNA) binding sites result from the fact that each TetR monomer binds to both DNA strands (5'-3' and 3'-5'). We will further investigate this in the subsequent section DNA Mutation.

Based on this, HDOCK computed the following output:

Figure 4: Result of the TetR-DNA docking simulation using HDOCK. The video was generated with PyMol.

Results

Figure 5: Transcription factor mutation alters DNA-binding affinity.

Figure 5 shows that most conservative mutations of the binding residues decrease the binding affinity. In case experimental results show that our quorum sensing based feedback loop gets activated at a low cell concentration not unique to solid tumors, we could substitute for instance the proline at position 39 with an alanine and therefore decrease the constant \(k_{T_{di}}\).

Only two conservative mutations (Ser at position 40 and Tyr at position 43) lead to an increased binding affinity. This indicates that the TetR system is already well optimized. We further investigated whether the binding affinity could be increased by introducing both mutations simultaneously to the TetR dimer. However, a TetR dimer with a Ser at position 40 and Tyr at position 43 in both chains binding to the TRE DNA motif results in an absolute HDOCK docking score of 203.73. This suggests that both mutations influence each other and therefore, when simultaneously introduced to the TF, slightly decrease the binding affinity.

In addition to the mutagenesis of our TF, we investigated the effect of the base composition of the DNA binding motif on the TF-DNA binding affinity. Figure 6 depicts the DNA binding sequence with the binding motif centered within the sequence, as described before. Both TetR monomers A and A' bind to both DNA strands. Furthermore, the DNA binding motif contains a sequence X , one guanine, and the reverse complement of X .

Figure 6: TRE DNA binding motif, based on Figure 3A in Ramos et al., 2005 [4].

Mutating the cytosine at position 9 in the 5' to 3' strand to an adenine therefore implies that the guanine at position 21 in the same strand is mutated to a thymine. Consequently, the reverse strand now also contains two base substitutions.

We investigated the effect of mutated bases at positions 9 to 14 in the 5' to 3' strand. As described, mutating one of these bases results in a total of four substituted bases. For each position, we substituted the respective base with the three other bases once. This results in 3 bases * 6 positions = 18 DNA strands. Subsequently, each DNA strand was docked to the baseline TetR dimer using HDOCK.

Results

Figure 7: TF-DNA binding affinity altered by base mutations introduced to the TRE DNA binding motif.

Figure 7 depicts the effect point mutations in the TRE DNA binding motif have on the TF-DNA binding affinity. At each position in the 5' to 3' DNA strand (x-axis), the original baseline expression is denoted as zero. Based on that, the other scores show the deviation from the baseline docking score (205.63, see Figure 5). For instance, the substitution of the cytosine at position 9 with an adenine results in an absolute docking score of 169.45 and \(169.45 - 205.63 \approx -36\). Therefore, a negative deviation in Figure 7 implies that the respective mutation decreased the TF-DNA binding affinity.

As already seen when investigating the effect of a TetR mutagenesis on the TF-DNA binding affinity, the results in Figure 7 suggest that the TetR-TER binding already is a well optimized system, because only one investigated mutation has the capability to increase the binding affinity: Substituting the adenine at position 14 with a cytosine results in a higher HDOCK docking score than the baseline.

Conclusions

Importantly, we are not exclusively aiming to increase the TF-DNA binding affinity. In order to fine-tune our feedback loop activation threshold, it is just as important to identify possible mutations that maintain the binding capability, but make it less energetically preferable. The performed investigation of effects of TetR as well as DNA mutagenesis, enables us to precisely optimize our system.

The following section describes an analysis of the docking of an extracellular ligand to two nanobodies (MESA) on the cell surface of our engineered CAR T-cells (see Figure 1, step 1). The rate constant used in our mathematical model to describe this process is \(K_{ML_e}\).

During the course of our project, we investigated two nanobodies: the GFP nanobody (GFPnb) and the mCherry nanobody (mCherryNB). Additionally, two ligands were examined: GFP-mCherry, with one GFP and one mCherry domain, as well as 2xmCherry, consisting of two mCherry domains.

As the structures of our constructs are unknown, we used AlphaFold [13] to predict a 3D model of both ligands as well as both nanobodies. Subsequently, we used HDOCK to simulate the docking between a nanobody dimer as a receptor and our ligand, consisting of two domains. For this, we identified the binding positions in the nanobodies by aligning our predicted protein structures with reference structures from PDB (3K1K for GFP; 6IR2 for mCherry).

We then performed three docking simulations, each with a different ligand - receptor combination.

Figure 10: HDOCK result of the docking simulation of the ligand 2xmCherry and two mCherry nanobodies as receptor. The video was created with PyMol.

Figure 10 displays one result of our docking simulation. The following table lists all performed simulations with their respective absolute HDOCK docking score.

The docking results correspond to the expected behavior. A mCherryNB dimer docked to the ligand 2xmCherry results in the highest docking score. Changing the ligand or one nanobody to GFB or GFPnb, respectively, decreases the absolute docking score. This implies a decreased docking affinity and consequently a decreased rate constant \(K_{ML_e}\).