Modeling: Functional optimization of the anticancer drug production
The targeted production of the anticancer drug 5-fluorouracil (5-FU) is the result of cascade of biomolecular reactions that begin with the anchoring of liposomes to cancer cells. As a trigger, the presence of the oncometabolite 2-hydroxyglutarate (2-HG) relieves gene repression exerted by DhdR [1]. In addition, activation of the T7 RNA polymerase leads to the transcription and subsequent translation of the thymidine phosphorylase gene [2]. The last step consists in the conversion of the prodrug Tegafur by the thymidine phosphorylase and the molecular diffusion of the product 5-FU toward the tumor.
To predict the dynamics of 5-FU production and support the experimental design, we built a kinetic model that describes the entire biomolecular network compartmentalized in a liposome. This model was used to demonstrate the feasibility of our project, predict the system’s behavior under various initial conditions, and design experiments. We carried out model-driven experiments and used our measurements to refine the model. Finally, we performed a sensitivity analysis to better understand how each species impacts 5-FU production, and we used this information to design an optimized version of our liposome with enhanced 5-FU production.
Model structure
Figure 1: Biomolecular network of the liposome represented following the Systems biology graphical notation (SBGN) [3]
The complete system was divided into three modules: (1) the cell biosensing module (in green) which encompasses the detection of the oncometabolite and the lifting of transcriptional repression, (2) the gene expression module (in blue) corresponding to the production of the enzyme thymidine phosphorylase from its gene, and (3) the anticancer drug production module that accounts for the enzymatic synthesis of 5-FU from the prodrug Tegafur (in orange). In this model, we consider that the split T7 RNAP has been refunctionalized following the anchoring of our liposome to a cancerous cell. We therefore assume that the polymerase is constitutively active.
Biosensing module
The IDH1 mutation present in many cancers leads to the accumulation of the oncometabolite 2-HG in the tumor's vicinity. This oncometabolite is present at relatively high concentrations in the tumoral environment and can be used as a disease biomarker to trigger a localized response.
We have thus designed a 2-HG biosensor based on a specific transcriptional regulator, DhdR. The DhdR protein may dimerize and interact with DNA through a specific sequence dhdO, thereby inhibiting gene expression. Repression can be lifted in the presence of 2-HG which binds to the dimer and induces its dissociation from DNA. Therefore, when the liposome enters the tumoral environment, 2-HG diffuses inside the lumen activating the expression of the drug-producing enzyme that was otherwise repressed by DhdR.
The mathematical representation of this 2-HG biosensor was inspired from iGEM team Duke 2021 [4]. The rates of all the reactions described below are modeled using a mass action law.
• Diffusion of DhdR into the liposome:
First, we used Fick's law to model 2-HG diffusion into the liposome:
- \(r_{2HG} = D_{2HG} ( 2HG_{out} - 2HG_{in} )\) (1)
Where D2HG is the permeability coefficient accross the liposomal membrane, and 2-HGin and 2-HGout are the concentrations of 2-HG inside and outside the liposome, respectively.
• Dimerisation of DhdR
- \(DhdR + DhdR\) ↔ \(D_2\) (2)
- \(r_{D2} = k_{D2} [DhdR] [DhdR] - k_{-D2} [D_2]\) (3)
Initially, the monomer DhdR can form a dimer, with kD2 and k-D2 the association and dissociation rate constants, respectively.
• Interaction of the DhdR dimer with dhdO
- \([dhdO - D_2]\) ↔ \([dhdO] + [D_2]\) (4)
- \(r_{DO} = k_{-D0} [dhdO - D_2] - k{DO} [dhdO] [D_2]\) (5)
The DhdR dimer D2 binds to the specific sequence dhdO with an association rate constant kDO and dissociation rate constant k-DO.
• Interaction of the DhdR dimer with 2-HG
2-HG can bind to free DhdR dimer, which prevents the repressor from interacting with DNA.
- \([D_2] + [2-HG]\) ↔ \([D_2 - 2HG]\) (6)
- \(r_{DH} = k_{DH} [D_2] [2HG] - k_{DH} [D_2 - 2HG]\) (7)
2-HG can also bind to the dimer that is already bound to DNA which results in the dissociation of the repressor-DNA complex.
- \([dhdO - D_2] + [2-HG]\) ↔ \([D_2 - 2HG] + [dhdO]\) (8)
- \(r_{DH} = k_{DH} [dhdO - D_2] [2HG]\) (9)
Here, kDH is the association rate constant of the DhdR dimer and 2-HG and k-DH the dissociation rate constant of the complex.
• Interaction of the T7 RNA polymerase with the promoter
- \([T7] + [dhdO]\) ↔ \([T7 - dhdO]\) (10)
- \(r_{DP} = k_{DP} [T7] [dhdO] - k_{DP} [T7 - dhdO]\) (11)
Once the split T7 RNA polymerase has been fully assembled after recognition of HER2, the enzyme can bind to the dhdO and activate the expression of tymp. Here, kDP is the association rate constant of dhdO and T7 RNAP and k-DP the dissociation rate constant.
Expression of the thymidine phosphorylase gene
The production of the Tegafur-activating enzyme Thymidine Phosphorylase (TYPH) is based on cell-free expression inside the liposome. In 2012, Stögbauer et al. [5] proposed a mathematical model to describe the expression of GFP in their cell free system. Therefore, we adapted this representation to the expression of TYPH in our PURE system by recalculating some parameters according to the length and sequence of our gene, as described below.
• Transcription
- \(\frac{d[mRNA]}{dt} = \frac{k_{ts}[TsR][dhdO - T7]}{Ks + [dhdO - T7]} - \delta_{mRNA}\) (12)
This equation represents the production of the tymp gene intro mRNA (using a law similar to the Michaelis-Menten equation) and the degradation or inactivation of mRNA using first-order kinetics. Ks is the Michaelis constant of transcription and kts the transcriptional rate constant.
As mentioned by iGEM team Delft 2021 [6], the transcription rate constant can be adapted to the length of the mRNA fragment depending on the composition of PURE system. The mRNA corresponding to the tymp gene is 1615 bp. According to Stögbauer et al., the transcription rate of the T7 polymerase (at 100 nM) is 2.2 NTP per second [5]. The constant kts can thus be calculated as:
- \(k_{ts} = \frac{NTP/s/polymerase}{NTP/mRNA}*[polymerase-in-PureFrex] = \frac{2.2}{1615}*0.1 = 1.35e^{-4}\) μ\(M.s^{-1}\)
In equation 12, [TsR] represents the amount of resources available in PURE system for transcription, i.e., of nucleotides. Since resources are finite, their consumption must be taken into account when modeling cessation of transcription due to resource depletion.
- \(\frac{[TsR]}{dt} = \frac{k_{cs}[TsR][dhdO - T7]}{Ks + [dhdO - T7]}\) (13)
Here, kcs represents the scaling factor for TsR. However, the value of kcs was determined for the GFP [5]. As we are expressing another protein with a different nucleotide sequence, we adapted the kcs value to our system based on the composition of PURE system, where all nucleotides are provided at the same concentration [7]. The limiting nucleotide of the gfp mRNA is adenosine with 245 adenosines per molecule. In our system, the limiting nucleotide for the tymp mRNA is guanine, with 472 guanines per mRNA. Hence, the resources are depleted 1.9 (\(\frac{472}{245}\)) times faster for the production of one tymp mRNA than for one gfp mRNA. We thus corrected kcs using the following formula:
- \(k_{cs} = 1.8e^{-4} * \frac{472}{245} = 3.39e^{3} s^{-1}\)
• Translation
The production of TYPH was modeled as:
- \(\frac{d[TYPH]}{dt} = \frac{k_{tl}[TlR][mRNA]}{Kl + [mRNA]}\) (14)
This equation represents the production of the enzyme TYPH from mRNA using the resources used for translation (TlR), i.e., the amino acids. Kl is the Michaelis constant of translation and ktl the translation rate constant.
As described above for kts, the value of the translation rate ktl was adapted to our protein sequence length which is 449 amino acids. According to Stögbauer et al. [5], the translation rate is 0.03 amino acid per second per ribosome, the latter being present at 2.5 μM in PUREfrex.
- \(k_{tl} = \frac{AA/s/ribosome}{AA/monomer}*[ribosomes] = \frac{0.03}{449}*2.5 = 1.67e^{-4}\) μ\(M.s^{-1}\)
Moreover, the translation resources TlR are used at a rate determined by the following equation:
- \(\frac{d[TlR]}{dt} = \frac{\delta_{TlR}[TlR]}{K_{TlR} + [TlR]}\) (15)
Here, KTlR is the Michaelis constant of TlR and δTlR the degradation rate. However, δTlR was only determined for GFP. Following the same reasoning as for NTP resources, we adapted δTlR to our system. While lycine is the limiting amino acid for GFP (22 Gly per protein), leucine is limiting for TYPH (70 Leu per protein). So, the consumption rate of the limiting amino acid was calculated as:
- \(\delta_{TlR} = 7.5e^{-5} * \frac{70}{22} = 2.39e^{-4} s^{-1}\)
Anticancer drug production and release
The third module describes the activation of the prodrug (Tegafur) into 5-FU. The enzymatic activity of the thymidine phosphorylase follows a Michaelis-Menten kinetic [8] and has already been characterized with Tegafur as the substrate.
- \(r_{5FU_{in}} = \frac{V_{max}[Tegafur]}{Km + [Tegafur]}\) (16)
The drug 5-FU produced in the vesicle limen will then diffuse outside the liposome:
- \(\frac{d[5FU_{out}]}{dt} = D_{5FU}( 5FU_{in} - 5FU_{out})\) (17)
Where D5-FU is the permeability coefficient across the liposomal membrane.
• No degradation of DNA or proteins (RNAP, TYPH, etc.).
• The thymidine phosphorylase follows a Michaelis-Menten kinetic.
• No evaporation.
• Diffusion is not limited inside the liposome, which is assumed to be homogeneous.
• The association rate of 2-HG with unbounded D2 is similar to the rate of 2-HG with bounded D2.
• 5-FU and 2-HG were considered as small polar molecules to approximate the corresponding membrane permeability coefficient.
Implementation
Overall, the model contains 16 species, 13 reactions, and a total of 25 parameters. The system of ordinary differential equations described above was implemented in COPASI (v4.38), and this model was used to simulate the temporal dynamics of all species. You can find this initial model in the Gitlab folder: Pre-calibration.
Model species and parameters
The list of all species and parameters of the model are given in this section.
Table 1 lists all species present in the liposome and their initial concentration.
Table 1: Species included in the model, with their initial concentration.
All model parameters are listed in Table 2 with the corresponding sources. We reused modeling works from previous iGEM projects: the dissociation and association rate constants related to DhdR and 2-HG were taken from iGEM Team Duke 2021 [4], and parameters related to gene expression were taken from iGEM Team Delft 2021 [6].
Table 2: Model parameters and values.
Results
We used this model to predict the dynamics of each species and analyse the drug production and its diffusion outside the liposome. The simulations obtained for a period of 6 hours are shown in Figure 2.
Figure 2: Evolution of the concentration of 2-HG in the liposome ([2-HGin], in orange), of the enzyme TYPH ([TYPH], in brown), of the 5-FU outside the liposome ([5-FUout], in red) and the 2-HG complexed to the dimer D2 ([2HG-D2], in green). Reaction time was 6 hours.
The model predicts that 2-HG diffuses inside the liposome and forms a complex with D2, hence decreasing the concentration of DhdR and lifting transcriptional inhibition. Consequently, the thymidine phosphorylase is produced rapidly, with a maximal concentration of 1.5 nM reached after 1h. Finally, Tegafur is converted to 5-FU, which diffuses outside the liposome.
The responsiveness of our anticancer liposome to a tumoral environment is a key design feature for limiting toxic side-effects. We thus used the model to evaluate whether 5-FU was produced at higher concentration in the presence of 2-HG (synthesized by cancer cells) compared to its absence. In absence of 2-HG, we expect the expression of TYPH to be repressed and, therefore, only a small amount of 5-FU would be produced. Consistently, simulations shown in Fig. 3 predicted that 1.5 μM of DhdR repressed gene expression, resulting in a ~2-fold lower concentration of 5-FU outside the liposome. Adding 100 μM of 2-HG, corresponding to typical concentrations in surrounding tumors [9], derepressed the system and boosted 5-FU production. The model therefore met the behavior expected for our liposomes and confirmed that our targeted drug production strategy could help tackle adverse effects of current chemotherapies.
Figure 3: Concentration of 5-FU after 6 hours in presence or absence of DhdR and 2-HG (concentrations are as indicated), with an initial Tegafur concentration of 1000 μM.
Guiding experiments and validating the model
We used the model to determine the optimal concentration of DhdR to be encapsulated in our liposome. To that aim, we built a reduced model containing only the biosensing and gene expression modules, along with a gene encoding the sfGFP protein that was used in our experiments as a reporter. You can find the reduced model in the Gitlab Folder: Experimental design.
As a first step, we predicted the sfGFP expression levels for different DhdR concentrations (between 0 and 1.5 µM).
Model predictions (line in Figure 4) indicate that, as expected, increasing DhdR concentration progressively reduces the amount of sfGFP produced.
Figure 4: Model predictions (in pink) qualitatively agree with experimentally measured fluorescence intensity values of sfGFP (in red) at different concentrations of DhdR.
Based on these encouraging predictions, we carried out the experiments with the same DhdR concentrations. A gradual decrease in the fluorescence intensity, i.e., in the expression of sfGFP, was experimentally verified (red bars in Figure 4). Importantly, the measurements are consistent with the model predictions both qualitatively and quantitatively confirming the parameter values and structure of the model. We decided to choose 1.5 μM as the optimal DhdR concentration for the next experiments as it efficiently repressed the expression.
As a second step, we wanted to validate the transcriptional activation by the 2-HG inducer. We thus predicted the expression level of sfGFP for a fixed concentration of DhdR (1.5 µM) and different concentrations of 2-HG (between 0 and 1 mM), and we carried out the corresponding experiments. Model predictions are compared to experimental results in Figure 5.
Figure 5: Model predictions (in pink) were compared to experimental fluorescence intensity values of sfGFP (in red) for a fixed concentration of DhdR (1.5 µM) and varying concentrations of 2-HG.
Experimental results confirm that increasing 2-HG concentration progressively lifts the repression caused by DhdR. However, the model poorly predicted the gradual increase of sfGFP levels measured experimentally, as it was less sensitive to lower concentrations of 2HG.
Calibration of the biosensing module using experimental data
We thus decided to use our measurements to redine our model. We estimated the values of different parameters of the biosensing module: kD2, k-D2, kDO, k-DO, kDH and k-DH by fitting our experimental data. For that, we used the Parameter Estimation task of COPASI (v4.38). Briefly, in this process, COPASI simulates the protein expression profiles for all the experimental conditions tested and estimates the parameter values that best fit the data (using the Particle Swarm algorithm as optimization method, with an iteration limit of 1000 and a swarm size of 40). You can find the experimental data and the calibrated model in the Gitlab Folder: Calibration. The best newly estimated values for the different model parameters are given in Table 3.
Table 3: Best newly estimated values for the different model parameters
In comparison to the predictions of the ab initio model constructed from the literature, the simulations with the refined parameter value were in quantitative agreement with our measurements (Fig. 6 and 7).
Figure 6: Comparison of the predictive capabilities of the kinetic model before (in pink) and after calibration (in green) using experimental data as inputs (in red). Concentrations of DhdR were as indicated and 2-HG was not present.
Importantly, the refined model with the calibrated biosensor module was also able to accurately simulate protein expression levels in the presence of 2-HG, corroborating its improved predictive power (Fig. 7).
Figure 7: Comparison of GFP expression levels predicted by the model before (in pink) and after calibration (in green) with experimental data (in red). Concentration of DhdR was fixed to 1.5 µM and different concentrations of 2-HG were used as indicated.
We then evaluated 5-FU production after 24 hours in response to different concentrations of 2-HG (Fig. 8) by comparing the simulation results before and after model calibration. You can find the model used in the Gitlab folder: Post-calibration.
Figure 8: Sensibility of the calibrated model (in red) to the concentration of 2-HG compared to the model before calibration (in pink). Reaction time was 24 hours.
Model simulations showed that the production of 5-FU after 24 hours was consistently lower with the updated model. Moreover, the yield of 5-FU production was very low, with a concentration of only 0.33 μM when 2-HG was present at concentrations encountered in the tumoral environment.
Enrico Mastrobattista, professor specialized in the preparation and characterization of Nanoparticles for Nucleic Acid and Protein Delivery into Cells, explained that liposomes are stable for 24 hours in a human body. Knowing that the IC50 of 5-FU, i.e., the concentration needed to inhibit a given biological process to half of the maximum, is approximately 46 μM [17], IC50 needs to be reached within those 24 hours.
Modeling results indicated that this concentration could not be reached within a reasonable time frame. We thus had to find a way to make our synthetic liposomes fully functional.
Based on the above modeling results, we wanted to determine how to optimize our system to reach IC50 within the 24 hour time frame, to ultimately propose an improved version of our liposome. Modifying kinetic properties of proteins (TYPH, T7 RNA polymerase, etc) is a challenging task that requires time-consuming enzyme engineering (e.g., there is no simple way to boost the enzyme's kcat). We therefore decided to explore a simpler and cheaper way to boost 5-FU production by focusing on the composition of our liposomes, hence opening a window towards features that we can easily control.
We applied a sensitivity analysis on the model to understand how the concentration of each species impacts 5-FU production. A sensitivity analysis provides a quantitative understanding on how a given parameter or concentration determines a given function of a system [18]. In our project, we define the function of interest as the time necessary to reach the IC50 of 5-FU outside the liposome (denoted 𝜏).
For each species, we calculated a control coefficient (C) that informs us on how its concentration (p) impacts 𝜏, using the following formula:
- \(C = \frac{(\tau(p)-\tau(1+\delta)p)/(\tau(p))}{(p-(1+\delta)p)/(p)}\) (18)
Each coefficient C quantifies the relative change in the time 𝜏 to reach IC50 in response to a relative change δ of the concentration p. If C is negative, an increase in p reduces 𝜏. On the contrary, if C is positive, an increase in p increases 𝜏. In our case, we want to decrease the time necessary to reach the IC50 for improving the cancer treatment. Therefore, we used the “Sensitivity analysis” task of COPASI to calculate the control coefficients of each species on 𝜏. Calculation results are shown in Table 4.
Table 4: Sensitivity analysis of the initial concentrations of each species on the time necessary to reach IC50. Green: Negative control coefficient, 𝜏 decreases when the concentration is increased. Red: Positive control coefficient, 𝜏 increases when the concentration is increased. White: The concentration does not impact 𝜏.
The results (Table 4) allowed us to identify the species that control 5-FU production, thereby guiding our liposome engineering strategy in a rational way. For example, the lowest coefficient is related to TlR, i.e., to the concentration of translational resources (amino acids) used for enzyme production. By increasing the concentration of specific amino acids in the liposome, especially leucine (the limiting amino acid), we expect a higher yield of 5-FU. We decided to quadruple the amount of resources in the liposome, from 1 μM to 4 μM, which is practically straightforward. In addition, since DNA, T7 RNA polymerase and transcriptional resources (TsR) also exert a significant negative control, we decided to double their concentration. Finally, we also sought to decrease the concentration of DhdR (from 1.5 μM to 1 μM) since this protein exerts a positive control. The proposed changes are summarized in Table 5.
Table 5: Initial and optimized concentrations based on sensitivity analysis results.
We then used the model to simulate the kinetics of 5-FU production with the refined concentrations (Fig 9). You can find the model in the Gitlab folder: Optimization.
Figure 9: Time course of 5-FU concentrations for the improved liposome designed based on the sensitivity analysis (red). The IC50 is appended in green. Reaction time was 24 hours.
This simulation indicates that, with the proposed modifications, the IC50 value of 5-FU can be reached in only 12 hours, which is compatible with the liposome stability in the body (about 24 h). This analysis is key to guide the design of new experiments with an improved composition of our therapeutic liposome.
Overall, the model showed that our synthetic liposomes can produce 5-FU at different concentrations depending on the presence of the oncometabolite 2-HG in the surrounding environment, hence demonstrating the feasibility of our project for targeted chemotherapy.
Based on these encouraging results, we then used the model to characterize in detail the biosensor module which is at the heart of our project. We designed and conducted model-driven experiments that quantitatively validated the strength of inhibition by DhdR as predicted by the model. In turn, we used these experimental data to improve our model, in particular transcription induction by 2-HG.
This improved model indicated that 5-FU production might not be sufficient yet to kill cancer cells. We thus carried out a sensitivity analysis which provided critical information that allowed us to design a new version of our liposomes. Specifically, we propose to adjust the concentrations of amino acids, DhdR, T7 RNA polymerase, and DNA, leading to a yield of 5-FU that reaches lethal concentrations in about 12 hours.
