Mathematical modeling plays a crucial role in understanding the project keeping precision l'esprit. We utilised modelling to gain an understanding of portions of the project we could not experimentally characterise. We aimed to build predictions that we could cross-validate with wetlab data. For Duonco, we segregated the entire project into its constituent steps to look at the processes we could model.
Figure 1: Project Plan
We began by identifying seven areas of the project we could work on. From our interactions with experts, we realistically narrowed down to five. These avenues have been described below.
One of the main questions we wanted to answer was: how could we predict the rate of production of our heterologous proteins inside the bacteria? This is important to predict density of the affibody and THP on Affi-OMVs and THP-OMVs respectively. To answer this question, we initially made a diagram of all the genes and proteins relevant to expression of the recombinant proteins from the plasmid. Our gene of interest is inserted into the E.coli expression vector pGEX-4T-1. Protein expression is under the control of the IPTG inducible tac promoter.
Figure 2. Gene circuit of our project
We wished to predict protein output (the dependent variable) and a function of IPTG concentration and time (independent variables). We applied the following assumptions while building our model:
Equations for the genomic copy are as given below:
(1.a) Constitutive transcription of LacI is modelled as with zero-order kinetics
\[ \phi \overset{k_{SMR} }\to M_{R}\]
\[ \frac{dM_{R} }{dt} = k_{sMR} - \lambda_{M_{R} } M_{R}\]
(1.b) Translation of LacI mRNA to the repressor protein R, follows first order kinetics
\[ M_{R} \overset{k_{SR} }\to M_{R} + R\]
(1.c) Dimerization of R follows second order kinetics
\[ 2R\underset{k_{-2R} }{\stackrel{k_{2R} }{\rightleftharpoons} }R_{2}\]
\[ \frac{dR}{dt} = k_{sR}M_{R} + 2k_{-2R} R_{2} - 2 k_{2R}{R}^{2} - \lambda_{R} R\]
Since the concentration of the operator is constant, it is equal to that of the combination of free operator and repressor bound operator
\[ O_{T} = O + R_{2}O\]
\[ \frac{dR_{2} }{dt} = k_{2R}{R}^2 - k_{-2R} R_{2} - k_{r} R_{2} O + k_{-r} (O_{T} - O) - k_{dr1} R_{2} {I}^{2} + k_{-dr1} [I_{2}R_{2}] - \lambda_{R_{2} } R_{2}\]
(2.a) The dimer \(R_{2}\) binds with the operator O in a 1-1 manner
\[ R_{2} + O \underset{k_{-r} }{\stackrel{k_{r} }{\rightleftharpoons} } R_{2}O\]
(3.a) 2 IPTG (I) molecules bind with 1 dimer \(R_{2}\)
\[ 2I + R_{2} \underset{k_{-dr1} }{\stackrel{k_{dr1} }{\rightleftharpoons} } I_{2}R_{2}\]
(3.b) IPTG can also bind to the repressor-operator\(R_{2}O\) complex such that 2 IPTG molecules bind to
1\(R_{2}O\)
\[ 2I + R_{2}O \underset{k_{-dr2} }{\stackrel{k_{dr2} }{\rightleftharpoons} } I_{2}R_{2} + O\]
\[ \frac{dO}{dt} = -k_{r} R_{2} O + k_{-r} (O_{T} - O) + k_{dr2} (O_{T} - O) {I}^{2} - k_{dr1} R_{2} {I}^{2}\]
(4.a) Transcription of the genomic LacY gene under the lac promoter is dependent on the concentration of free operator
\[ O \overset{k_{S1MY} }\to O + mRNA_{Y}\]
\[ \frac{dM_{Y} }{dt} = k_{s0MY} (O_{T} - O) + k_{s1MY} O - \lambda_{MY} M_{Y}\]
(4.b) Leaky transcription of LacY can be modelled as dependent on concentration of bound receptor \(R_{2}O\)
\[ R_{2}O \overset{k_{SOMY} }\to R_{2}O + mRNA_{Y}\]
(4.c) Translation of LacY protein(Y) follows first order kinetics
\[ M_{Y}\overset{k_{SY} }\to M_{Y}+Y\]
(5.a) The facilitated tranportation of IPTG by Y from the extracellular region where its concentration is \(I_{ex}\)
to the
intracellular region where its concentration is I. This reaction is modelled with Michaelis Menten equations
\[ Y + I_{ex} \underset{k_{-p} }{\stackrel{k_{p} }{\rightleftharpoons} }YT_{ex} \stackrel{kft}\to Y + I\]
\[ \frac{d YI_{ex} }{dt} = -(k_{ft} + k_{-p}) YI_{ex} + k_{p} Y I_{ex} - \lambda_{YI_{ex} YI_{ex} }\]
\[ \frac{dY}{dt} = k_{sY} M_{Y} + (k_{ft} + k_{-p}) YI_{ex} - k_{p} Y I_{ex} - \lambda_{Y} Y\]
(5.b) IPTG diffuses across the membrane
\[ I_{ex} \underset{k_{t} }{\stackrel{k_{t} }{\rightleftharpoons} }I_{in}\]
\[ \frac{dI_{in} }{dt} = -2 k_{dr1} R_{2} {I}^{2} + 2 k_{-dr1} [I_{2}R_{2}] - 2 k_{dr2} (O_{T} - O) {I}^{2} + 2 k_{-dr2} O [I_{2}R_{2}] + k_{ft} YI_{ex} + k_{t} (I_{ex} - I) + 2 \lambda_{[I_{2}R_{2}]} [I_{2}R_{2}] + \lambda_{YI_{ex} } Y_{I_{ex} }\]
(6.a) The degradation equations for various molecular species are denoted below:
\[ M_{R}\overset{\lambda_{M_{R} } }\to\phi\]
\[ M_{Y}\overset{\lambda_{M_{Y} } }\to\phi\]
\[ R\overset{\lambda_{R} }\to\phi\]
\[ R_{2}\overset{\lambda_{R2} }\to\phi\]
\[ Y\overset{\lambda_{Y} }\to\phi\]
\[ YI_{ex}\overset{\lambda YI_{ex} }\to I\]
\[ I_{2}R_{2}\overset{\lambda_{I_{2}R_{2} } }\to I\]
\[ mRNA_{Y}\overset{\lambda M_{Y} }\to \phi\]
\[ C\overset{\lambda c}\to \phi\]
(7.a) Transcription of LacI under the LacI promoter in the plasmid
\[ \phi \overset{k_{SMR} }\to M_{R}\]
\[ \frac{dM_{Rp} }{dt} = k_{sMR} - \lambda_{M_{R} } M_{Rp}\]
(7.b) Translation of LacI mRNA to the repressor protein, denoted Rp for the plasmid
\[ M_{Rp} \overset{k_{SR} }\to M_{Rp} + Rp\]
(7.c) Dimerization of Rp for the plasmid
\[ 2Rp\underset{k_{-2R} }{\stackrel{k_{2R} }{\rightleftharpoons} }R_{2p}\]
\[ \frac{dRp}{dt} = k_{sR}M_{Rp} + 2k_{-2R} R_{2p} - 2 k_{2R}{Rp}^{2} - \lambda_{R} Rp\]
Here too, the concentration of the operator is constant, it is equal to that of the combination of
free operator and repressor bound operator
\[ O_{Tp} = Op + R_{2p}Op\]
\[ \frac{dR_{2p} }{dt} = k_{2R}{Rp}^2 - k_{-2Rp} R_{2p} - k_{r} R_{2p} Op + k_{-r} (O_{Tp} - Op) - k_{dr1} R_{2p} {Ip}^{2} + k_{-dr1} [I_{2}R_{2}p] - \lambda_{R_{2} } R_{2p}\]
(8.a) The dimer \(R_{2p}\) binds with the operator Op of plasmid in a 1-1 manner
\[ R_{2p} + Op \underset{k_{-r} }{\stackrel{k_{r} }{\rightleftharpoons} } R_{2p}Op\]
(9.a) IPTG(denoted I) binds with the dimer R_{2p} such that 2 IPTG molecules bind with 1 dimer \(R_{2}\)
\[ 2I + R_{2p} \underset{k_{-dr1} }{\stackrel{k_{dr1} }{\rightleftharpoons} } I_{2}R_{2}p\]
(9.b) IPTG can also bind to the repressor-operator\(R_{2}Op\) complex such that 2 IPTG molecules bind to
\(R_{2}Op\)
\[ 2Ip + R_{2}Op \underset{k_{-dr2} }{\stackrel{k_{dr2} }{\rightleftharpoons} } I_{2}R_{2}p + Op\]
\[ \frac{dOp}{dt} = -k_{r} R_{2p} Op + k_{-r} (O_{Tp} - Op) + k_{dr2} (O_{Tp} - Op) {I}^{2} - k_{dr1} R_{2p} {I}^{2}\]
(10.a) Transcription of ClyA(denoted c) under the pTac promoter is dependent on concentration of free operator
and the rate constant is 11 times \(k_{S1MY}\)
\[ O \overset{k_{S1MYp} }\to O + M_{Yp}\]
\[ \frac{dM_{Yp} }{dt} = k_{SIMPp} O + K_{leakyp}\]
(10.b) Leaky transcription of ClyA
\[ R_{2p}Op \overset{k_{leakyp} }\to R_{2p}Op + M_{Yp}\]
(11.a) Translation of \(M_{Yp}\) to the protein ClyA
\[ M_{Yp} \overset{KM_{Yp} }\to M_{Yp} + C\]
\[ \frac{dC}{dt} = k_{M_{Yp} } M_{Yp}\]
| Symbol | Species denoted |
|---|---|
| MR | LacI repressor mRNA |
| R | LacI repressor monomer |
| R2 | LacI repressor dimer |
| O | lacO operator |
| R2O | Repressor-operator complex |
| I | Intracellular inducer IPTG |
| Iex | Extracellular inducer IPTG |
| I2R2 | Repressor-inducer complex |
| MY | LacY permease mRNA |
| Y | LacY permease |
| YIex | Permease-inducer complex |
| C | Cly-A Protein complex |
| \(\phi\) | Generic source or sink |
The constants and their values used are given in a table here.
| Symbol | Value | Units | Description |
|---|---|---|---|
| VE. coli | 8 × 10−16 | L | E.coli Volume |
| OT | 1 | (copy number) | Copy number |
| [OT] | 2.08 nM | nM | Concentration of operator |
| ksMR | 0.23 | nM·min−1 | lacI transcription rate |
| ksR | 15 | min−1 | LacI monomer translation rate constant |
| k2R | 50 | nM−1·min−1 | LacI dimerisation rate constant |
| k−2R | 10−3 | min−1 | LacI dimer dissociation rate constant |
| kr | 960 | nM−1·min−1 | association rate constant for repression |
| k−r | 2.4 | min−1 | dissociation rate constant for repression |
| kdr1 | 3 × 10−7 | nM−2·min−1 | association rate constant for 1st derepression mechanism |
| k−dr1 | 12 | min−1 | dissociation rate constant for 1st derepression mechanism |
| kdr2 | 3·10−7 | nM−2·min−1 | association rate constant for 2nd derepression mechanism |
| k−dr2 | 4.8·103 | nM−1·min−1 | dissociation rate constant for 2nd derepression mechanism |
| ks1MY | 0.5 | min−1 | lacY transcription rate constant |
| ks0MY | 0.01 | min−1 | leak lacY transcription rate constant |
| ksY | 30 | min−1 | lacY translation rate constant |
| kp | 0.12 | nM−1·min−1 | LacY-IPTGex association rate constant |
| k−p | 0.1 | min−1 | LacY-IPTGex dissociation rate constant |
| kft | 6.104 | min−1 | IPTG-facilitated transport constant |
| kt | 0.92 | min−1 | IPTG passive diffusion constant |
| kMc | 30 | min−1 | Rate of translation of ClyA |
| kS1Myp | 0.5 | min−1 | Rate of transcription of ClyA |
| λMR | 0.462 | min−1 | lacI mRNA degradation constant |
| λMY | 0.462 | min−1 | lacY mRNA degradation constant |
| λR | 0.2 | min−1 | repressor monomer degradation constant |
| λR2 | 0.2 | min−1 | LacY degradation constant |
| λY | 0.2 | min−1 | LacY-inducer degradation constant |
| λYIex | 0.2 | min−1 | LacY-inducer degradation constant |
| λI2R2 | 0.2 | min−1 | repressor-inducer degradation constant |
| λIin | Negligible | min−1 | Degradation constant of IPTG inside cell |
Figure 3: Protein concentration vs IPTG concentration
We obtained a plot of indicating the relationship between the concentration of the inducer IPTG and the output protein. We found that the curve saturates ~1 mM IPTG, and utilised this value while inducing our own bacterial cultures. We found sufficient protein production following this.
One of the processes we aimed to modeling OMV internalisation. OMV internalisation is defined as the process of OMVs entering mammalian cells by receptor mediated endocytosis. In this case, the receptors are HER2 and CX3CR1. Building a model which can recapitulate the essential features of the process can help us predict the concentrations of the enzyme and prodrug to load. This model can also be used to gain estimates on what percentage of administered OMVs internalise, and what the the maximum dosage of OMVs should be.
The following equations can be written to describe the process
| Symbol | Species denoted |
|---|---|
| \(C\) | Concentration of Vesicles |
| \(C_{b}\) | Concentration of bounded vescicles |
| \(C_{bs}\) | Concentration of free binding sites |
| \(C_{i}\) | Concentration of internalised vesicles |
Solving the above system of equations numerically gives us the plots below:
Figure 4: Internalisation of OMVs at low concentrations follows a linear pattern
Electroporation is the method of using high voltage electric pulses to increase cell permeability facilitating delivery of drug and genetic material into the cell. In our project, this technique plays a crucial role in loading the prodrug into the OMVs. The challenge was to optimise the size of the pore so as to enhance the rate of drug carried in each OMV without compromising the drug loading process or membrane integrity. Hence, a mathematical model was created to calculate the optimum pore size for effective cargo loading. This was done with the help of one-dimensional differential equations that took the pore size and carrying capacity as the variables.
We assume the OMVs to be spherical in nature, occupying a radius R. The electric field applied induces pores on the surface. Their exact size and radius distribution depends on the electric field being applied. The pore density N is governed by the well known Krassova and Filev equation, given below \[ \frac{d N}{d t} = \alpha\exp\left(\frac{V_{m} }{V_{e p} }\right)^{2}\left[1-\frac{N}{N_{0}\exp\left(q\left(\frac{V_{m} }{V_{e p} }\right)^{2}\right)}\right]\]
Where {V_{m} } is the transmembrane potential of the spherical cell in a uniform electric field defined as: \[ {V_{m} } = 1.5E \times r_{c} cos \psi\]
We isolate particulars of these equation using the below constants \[ K = \alpha{e}^{\frac{V_{m}^{2} }{V_{ep}^{2} } }\] \[ C = \frac{ {e}^{-q{\frac{V_{m}^{2} }{V_{ep}^{2} } } } }{N_{o} }\]
The equation then becomes:
\[ \frac{dN_{ep} }{dt} = K(1 - N_{ep}C)\]
Solving the above ordinary differential equation yeilds us:
\[ -{\frac{ln(1 - N_{ep}C)}{C} }{\rvert}_{N_{epo} }^{ {N_{ep} }^{t} }\]
\[ -KCt = ln{(\frac{1 - {N_{ep} }^{t}C}{1 - {N_{ep} }^{o}C})}\]
\[ e^{-kct} = \frac{1 - {N_{ep} }^{t}C}{1 - {N_{ep} }^{o}C}\]
\[ e^{-kct}{(1 - {N_{ep} }^{o}C)} = {1 - {N_{ep} }^{t}C}\]
The constants for our case are given below:
\[ \alpha = 1\times10^{9}\]
\[ V_{m} = 0.08\]
\[ V_{ep} = 0.258\]
\[ {(\frac{V_{m} }{V_{p} })} = 0.3100775194\]
\[ {(\frac{V_{m} }{V_{p} })}^{2} = 0.096148\]
\[ e^{ {(\frac{V_{m} }{V_{p} })}^{2} } = 1.10092\]
\[ e^{-{(\frac{V_{m} }{V_{p} })}^{2} } = 1.10092\]
\[ r_{m} = 0.8\times{10}^{-9}m\]
\[ r_{a} = 0.51\times{10}^{-9}m\]
\[ \frac{r_{m} }{r_{\alpha} } = 1.5686\]
\[ q = 2.46059\]
\[ K = \alpha{e}^{\frac{V_{m}^{2} }{V_{ep}^{2} } } = 1.10092\]
\[ C = \frac{ {e}^{-q{\frac{V_{m}^{2} }{V_{ep}^{2} } } } }{N_{o} }= {(\frac{0.789322}{N_{o} })}\]
We plug in these values for our original equation. It leads us to the plot shown below.
\[ {N_{ep} }^{t} = {\frac{1 - e^{-kct}{(1 - {N_{ep} }^{o}C)} }{C} }\]
Figure 6: Pore density evolution with time. As expected, the pore density saturates off to a particular value depending on the field strength(characterized by membrane voltage Vm). One also expects the pore density to exponentially decay off to zero as we remove the force field.
The prodrug is the precursor drug that is inactive but, when acted upon by the cognate enzyme, is activated. The damage to healthy cells is prevented by this inactive nature of prodrug, as the specific
combination of HER2 and CX3CR1 found on cancerous cells will make sure that the prodrug and enzyme are selectively internalised. The prodrug-enzyme kinetics are
important in understanding amount of prodrug and enzymes required to trigger cell death.
The prodrug in our project is 5-fluorouracil (5 FU) and the enzyme is cytosine deaminase. 5 FU is activated to 5-fluorocytosine by cytosine deaminase.
| Rate | \[ {k_{1} }^{app}\] | \[ {k_{3} }^{app}\] | \[ \frac{k_{1} \times k_{3} }{k_{2}+k_{3} }\] | \[ {k_{1} }^{app}\] |
|---|---|---|---|---|
| 1 | 0.316 | 0.0029 | 0.0012 | 0.0090 |
| 2 | 0.831 | 0.0023 | 0.0025 | 0.0018 |
| 3 | 0.728 | 0.0010 | 0.0010 | 0.0009 |
| 4 | 0.287 | 0.0022 | 0.0008 | 0.0010 |
| 5 | 0.267 | 0.0028 | 0.0010 | 0.0011 |
\[ k_{i} = \frac{ {k_{1} }^{app}\times{k_{3}^{app} } }{ {k_{2} }^{app} + {k_{3}^{app} } }\] \[ Rate = -\frac{d[A]}{dt} = k[A]\] \[ \frac{d[A]}{A} = -kdt\] \[ \int_{[A]_{o} }^{[A]}\frac{d[A]}{[A]} = -\int_{to}^{t}kdt\] \[ \int_{1}^{[A]}\frac{d[A]}{[A]} = -\int_{to}^{t}kdt\] \[ ln[A] - ln[A_{o}] = -kt\] \[ ln[A] = -kt + ln[A_{o}]\] \[ t_{1/2} = \frac{ln2}{k}\] \[ t_{1/2} = 551.26\; \pm 404.92s\]
We calculated the half-life of 5-FU to be 551.2 s
Ligand-receptor interactions are a major class of protein-protein interactions that regulate life
processes and helps maintain specificity of various proteins and their functions. The conformational
change that the ligand undergoes effectively acts as a lock and key mechanism which discriminates
between proteins and makes sure that only the specific ligand binds to the receptor site. In the context of our project, the ligand is the affibody or the THP and the receptors are HER2 and CX3CR1 respectively.
Understanding ligand-receptor interactions and quantifying their rates are crucial for designing the prodrug delivery mechanism and ensuring the selective targeting of cancerous cells with minimum
damage to healthy cells. The rate at which the ligand binds to the receptor, the mechanism of binding and the
probability of a ligand failing to bind give us important information in deciding the optimum dosage of the prodrug to be loaded. It also helps in finding the rate of reaction and the time taken for the
successful completion of the binding. Such studies can be done in-silico with Autodock Vina,
a popular molecular docking software.
We went on to conduct docking studies for our project. One of the main inputs we required was the structure file of our THP. Since we were unable to find the THP structure anywhere, a protein structure prediction was done using AlphaFold, an artificial intelligence program developed by DeepMind for the prediction of protein structure. The predicted structures are shown below.
Figure 6: Structures of THP
We required a ligand-receptor system to perform the docking and the following molecules coupled with the above molecule made the system complete.
Figure 7:Other Proteins: Affibody(Left), CX3CR1(Center), HER2(Right)
We attempted docking using Autodock Vina. However, after multiple attempts over multiple days, we were unable to proceed with the same. After troubleshooting the issue, we came to the conclusion that the molecules we used were too big for the software to handle, thus, our dream of a theoretical free energy estimate was lost. This problem was rectified using literature sources to find our k values.
