Team:Amazonas-Brazil/Model
Model
Disclaimer: If you are already familiar with modelling and its tools, skip to the Our model section.
Introduction
What is modelling?
Math is grounded in axioms based on logic, which is the most basic way one may achieve conclusions. These axioms construct a system that always needs to respect them. So, math is as deep as one can get, and have results from logic assumptions.
Bearing this in mind, it's only natural that this axiomatic logic can be used as a basis for constructing more and more complex things. But what is modelling?
Modelling describes a system through assumptions that can be interpreted as logic respecting the proposed axioms. This means that when a certain set of axioms is chosen, all the tools developed for it can also be used for the constructed model.
This line of reasoning is crucial for modern science because a good model can encapsulate every important aspect of a system through its assumptions, granting a way to get data without experiments! The first science to use it successfully was physics, developing ways to predict behaviors otherwise unknown, going beyond lab data from only math itself.
How to build a model?
As explained earlier, a model consists of simplifications of reality, in a way possible to describe it mathematically. So, first of all, we need to make assumptions to simplify and turn the system into something treatable with a certain mathematical framework enabling all of its tools.
This seems pretty abstract, but is actually very intuitive and can be better understood with an example. Let's say you want to build a model to describe winds in a certain place. Wind currents have a velocity, also a predominant direction and orientation... Very interesting way to think about it, because it looks just like vectors! Where the speed can be represented as the module of the vector. Not only that, but every point in space has a vector attached to it, meaning we have a vector field! (See Figure 1)
Figure 1: Vector field representing the speeds of wind above 40m in Southern California[1]. Available via license:
Creative Commons Attribution 4.0 International
After choosing the framework and making necessary assumptions to treat this problem, all of the vector field treatment tools are enabled for us. For example, there is an algebraic topology theorem, the hairy ball theorem, which states that in a vector field tangent to the surface of a sphere, there needs to exist at least one point where the vector is zero. Extending the vector field treatment to the whole earth, this suggests that at all times, wind (our vector field) must have at least one zero point, in this case, a zero-point represents a cyclone/anticyclone, like the ones observed in Figure 2. Which seems strange but at the same time is incredible! From pure math and assumptions, we come to a true observable conclusion: at all times there must exist at least one cyclone/anticyclone on the earth! Even though we can't know exactly where it exists, we know it does.
Figure 2: Gif showing how each of the behaviors of a cyclone/anticyclone and saddle that can occur in a zero wind place. Source:
http://chalkdustmagazine.com/blog/hairy-balls-cyclones-computer-graphics/
This simple example shows us how a model works and predicts behaviors, arriving at conclusions otherwise unknown, without experiments. Not only that, but models can help the experiments to be more precise and also expand the data beyond measurements, guaranteeing a certain behavior, the same way that experiments validate and enhance models, turning into a feedback loop, where experiments feed the model and vice-versa. Furthermore, models also turn works more scientific, showing that the experimental data can be expected, turning it more credible and reproducible, as we now know what aspects are essential for it to work.
Models in Synthetic Biology
In Synthetic Biology, experimentation can be pretty expensive, in an ideal case, we want to do the least experiments possible. This aspect is what makes models essential for this new scientific branch, cutting costs, predicting behaviors, helping characterization and expanding data overall. This can be noted from the works of the repressilator[2] and the Collins toggle switch[3], that heavily used mathematical modelling in their works, showing the expected performance.
Mathematical tools
This subsection aims to give a brief introduction to the math used in our model, which includes: logic gates and dynamical systems, also talking a bit about the construction of circuits with the help of Hill function, allowing the unfamiliar to keep up with the work developed here. Synthetic Biology draws lots of tools from electronics and engineering such as logic gates and dynamical systems which are used to describe circuits with varying complexity.
\(\boldsymbol{\cdot} \) Logic gates
In electronics, circuits can have an unending amount of complexity added to it, in a way that it can get quite convoluted and almost impossible to take information from it. Logic gates come to simplify the behaviors implementing a Boolean algebra, where 1 represents the presence of a certain component, and 0 the absence. The most common are the: NOT, AND, OR, NAND, NOR, XOR and XNOR gates, represented in Figure 3.
Figure 3: All logic gates representations possible, where these symbols are implemented to build a certain circuit. Source: https://physicsabout.com/logic-gates/.
For example, let's take a look at the NOT gate. The NOT gate means that when there is no input, we have an output, and where there is input, there is no output. For example, every time it snows you don't go to school, snow is the input and going to school is the output, and the relation between going to school and snowing can be seen as a NOT gate. Even though this is a very simple example, this tool is very important for the construction of complex circuits, as one might look at the funny symbol and know what to expect from it, even if we don't see the full machinations behind it.
| Snow | Going to school |
|---|---|
| 0 | 1 |
| 1 | 0 |
Table 1: Demonstration of how the NOT gate works utilizing Boolean algebra for the example given above.
\(\boldsymbol{\cdot} \) Dynamical Systems
The first thing to note before starting to talk about dynamical systems is that I'm assuming the reader already has an understanding of basic calculus and its more "physical" interpretations. From calculus, it's important the notion that a rate of change can be described as a derivative, as can be seen from the Newton quotient. This idea is the foundation for this whole branch of mathematics known as Dynamical Systems.
To understand this, you can think of how speed is the derivative of distance with relation to time meanwhile, acceleration is the derivative of speed \(\vec{a}=\frac{d\vec{v}}{dt}\). In this way, you see the speed as the rate of change of the position of a certain object, which is very intuitive. Such intuition is so important that every science and engineer should understand it, as it is the starting point for most mathematical formalizations.
The cases described above are very simple, but gives us an idea of how these equations behave, other rates can also be considered, adding complexity. Let's say that a raindrop is falling, but as it gains speed, it accumulates more water and its radius becomes bigger by a factor k and g is the gravity. So we can describe the system as:
$$ \frac {dv}{dt} = g$$
$$ \frac {dR}{dt} = kv$$
Even though it is not physically correct, it gives us a feeling of how these ideas work. This particular dynamical system can be solved analytically, but the majority actually can't. So why study something we can't solve? One of the most important tools for dynamical systems is the phase portrait or phase diagram. Such a tool allows us to have more qualitative behaviors, without an analytical solution, predicting if the solutions will oscillate, be attracted by a particular point or repelled, and many more interesting behaviors [4]. The way it works is that we implement the rates as vector in a vector field, where the axis are the analyzed quantities, as seen in Figure 4.
Figure 4: Phase portrait of a pendulum, where its axis are the value theta of the angle and the angular velocity, where the balls represent a specific initial condition and latter how it evolves with time. Source: https://gereshes.com/2019/03/04/an-introduction-to-phase-portraits/
Even though there are other ways to mathematically describe circuits created from Synthetic Biology, such as the thermodynamical or the stochastic approach, the Dynamical System way is for sure the most widely used, being essential for the Synthetic biologist to understand the basics.
Model construction
Now that we have an idea of how Dynamical Systems work, we can start building the mathematical description of simple circuits and the assumptions needed.
The first and more important assumption is that there are enough copies of the analyzed component, in a way that a continuous approximation is good. When dealing with Synthetic Biology, in general, there are enough copies to make this consideration, as the colonies tend to be big enough. Another important aspect that in general is used is that a rate instantly interferes with the others, such that there are no delays from one to another. More complex models may consider this, with discrete-time lags, but it will not be considered further in this work. Also, this type of models doesn't consider how the reactions occur, only how their rate change and affect other rates.
With these assumptions, a simple model can be built for a chemical reaction such that:
$$ E + S \rightarrow ES $$
This can be mathematically described as, where the brackets represent the concentration and k is the reaction constant:
$$ \frac{d[E]}{dt} = -k[E][S] $$
$$ \frac{d[S]}{dt} = -k[E][S] $$
$$ \frac{d[ES]}{dt} = k[E][S] $$
As more and more reactions occur, more terms appear in each equation, the positive ones are the production rates, while the negative ones are the consumption rates. These equations can be treated as previously with dynamical systems, allowing us to describe how each concentration of E, S and ES change with time.
Another important tool for Syn Biology is the Hill equation, observed in Figure 5, where it usually appears every time there is activation or repression of a certain concentration. For example, suppose a concentration of A regulated by L, the Hill equation is written as for activation:
$$ \frac{dA}{dt} = \frac{v_{max}[L]^ {n} }{ K + [L]^n} $$
And for repression:
$$ \frac{dA}{dt} = \frac{v_{max} }{ 1 + \frac{[L]^ {n}}{K} } $$
Figure 5: Representation of the Hill function for activation, where we can see its sigmoidal shape[5].
Where \(v_{max}\) is the maximal production rate, K is the Michaelis Menten constant, which tells us about the half-maximum production and n is the Hill coefficient, which is related to the steepness of the curve, or how fast it goes to its minimal or maximal value. Another interesting thing to notice is that for the first equation, as L grows, the rate goes to \(v_{max}\) in the limit. And for the second equation, the rate is \(v_{max}\) when [L] = 0 and goes to zero as [L] grows [6]. This is exactly the behavior wanted from a repression and activation function.
Our model
Our model is divided in 3 major parts:
The first one helped our design to choose the topology that minimized leakiness, as it is essential for our system to focus on tumor damage only.
The second part was a validation to go beyond experimental data, provinding us more insight into what we built with the help of a modelling as described in the other subsection, using deterministic equations.
The third and last part was to predict how our system would impact tumor size, allowing us to give some idea of how the treatment would work if done in real life.
Design influence
Logic gates are schemes widely used in electronics that simplify complex inner components with a symbol that summarizes it efficiently with the help of Boolean algebra. This makes it easier to take information, even if we are working with a convoluted system with many moving parts.
Our system aims to create a logic gate AND that expresses a target protein only when two inputs (related to tumor micro ambient) are present, making a scheme as bellow:
| Lactate | Hypoxia | Output |
|---|---|---|
| 0 | 0 | 0 |
| 1 | 0 | 0 |
| 0 | 1 | 0 |
| 1 | 1 | 1 |
Table 1: Table demonstrating the way we want our design to work, utilizing Boolean algebra, where 1 represents the presence of the component, while 0 represents its absence.
Even though we know that biological systems can't have perfect signals, meaning that even if the Boolean algebra states that there is no output, there is a residual production, due to thermodynamical fluctuations, the same way that there is a growth curve when our system is "ON". These factors are crucial to our design because a great quantity of output in healthy areas can cause significant damage to healthy tissues, the same way that low production of target protein in the tumor tissues indicates that the treatment isn't effective.
Figure 6: Graph showing how a perfect AND gate behaves, where there is production only when the system is totally ON. We also observe that when it is ON, there is a discontinuity, indicating that the derivative in the point is infinity, something that we aim with our design.
In our case it is impossible to exist discontinuities, like the one in Figure 6, because we are treating a continuous model, however, we wish to construct a design that maximizes the derivative in the "switch ON point". We also want to minimize the OFF residual production and maximize the ON production. Given the inputs, we built two probable designs, which will serve as placeholders to future discussions, and through a deterministic model we will choose the most efficient one for our specifications.
The first circuit follows the logic bellow, where the input promoters of Lactate and hypoxia, express hypothetical proteins A and B, which form a complex AB that activates the output promoter.
Figure 8: Representation of circuit 1.
The second circuit follows a logic based on repression, where the input promoters produce hypothetical proteins that both repress the other two promoters, which repress the output promoter. The input promoters repress the repressor of the output promoter.
Figure 7: Representation of circuit 2.
Now that we have the representation, a system of deterministic equations can be built which we are going to analyze, constructing it as described in the introduction page.
Equations circuit 1
$$\frac{d[A]}{dt} = y_1 + \frac{(Y_1 - y_1)[Lact]^{n_1}}{ K_1 + [Lact] ^ {n_1}} - \delta_1 [A] \tag{1.a}$$
$$\frac{d[B]}{dt} = y_2 + \frac{(Y_2 - y_2)K_2}{ K_2 + [O_2] ^ {n_2}} - \delta_2 [B] \tag{1.b}$$
$$\frac{d[AB]}{dt} = \alpha [A][B] - \delta_3 [AB] \tag{1.c}$$
$$\frac{d[GFP]}{dt} = y_4 + \frac{(Y_4 - y_4)[AB]^{n_4}}{ K_4 + [AB] ^ {n_4}} - \delta_4 [GFP] \tag{1.d}$$
Equations circuit 2
$$\frac{d[A]}{dt} = y_1 + \frac{(Y_1 - y_1)[Lact]^{n_1}}{ K_1 + [Lact] ^ {n_1}} - \delta_1 [A] \tag{2.a}$$
$$\frac{d[B]}{dt} = y_2 + \frac{(Y_2 - y_2)K_2}{ K_2 + [O_2] ^ {n_2}} - \delta_2 [B] \tag{2.b}$$
$$\frac{d[C]}{dt} = 2y_3 + \frac{(Y_3 - y_3)K_3}{ K_3 + [A] ^ {n_3}} + \frac{(Y_3 - y_3)K_4}{ K_4 + [B] ^ {n_4}} - \delta_3 [C] \tag{2.c}$$
$$\frac{d[GFP]}{dt} = y_4 + \frac{(Y_4 - y_4) }{ K + [C] ^ {n_5}} - \delta_4 [GFP] \tag{2.d}$$
Explanations and manipulations
Where \(y_i\) and \(Y_i\) represents the minimal and the maximal expression of the ith promoter, A and B are hypothetical proteins that form a complex to activate a promoter, \(n_i\) represents the Hill coefficient of the ith promoter and K is the Michaelis Menten constant and \(\delta\) is the degradation + dilution term. Most equations were modelled using Hill's equation for activation and repression, while equation 1.c was modelled as a law of mass action, ignoring the backward reaction, as explained in the introduction page.
The first thing we’ve noticed from equations a and b in both circuits is that they have the following form:
$$\frac{d[A]}{dt} = \beta - \delta [A] \tag{3}$$
Where \(\beta\) is a constant. This means they can be solved analytically wielding results:
$$ A = -\frac{e ^ {-\delta (t + c_1)}} {\delta} + \frac{\beta}{\delta} \tag{4}$$
This is very interesting because we will find that A and B are constants when t goes to infinity, both can be substituted in equation 2.c, which can then be solved using the same method, as everything is constant, then substituted in equation 2.d, which can be solved to find the expression of GFP as time is big enough, depending on Lactate and Oxygen.
Now, we want to compare the leakiness of both circuits, meaning that we first want to compare when the system is fully off, meaning that \([O_2] \rightarrow \infty\) and \([Lac] \rightarrow 0 \) and for simplicity sake, we will assume that all promoters have the same associated parameters. Solving the ODEs as explained wields results:
$$ A = \frac{y} {\delta} \tag{5}$$
$$ B = \frac{y} {\delta} \tag{6}$$
This result can be substituted in the following equations meaning that we will have the following leakiness for each circuit:
$$\eta _1 = \frac{1}{\delta} \left(y + \frac{ (Y - y) ( \frac{\alpha y^2}{\delta ^3 })^ {n} } {K + ( \frac {\alpha y^2} {\delta ^3} ) ^ {n} } \right)\tag{7}$$
$$ \eta _2 = \frac{1}{\delta} \left( y + \left( \frac{(Y-y)} { K + 2 \left(y + \frac{(Y-y)}{K + (\frac{y}{\delta}) ^{n}} \right) ^{n} } \right) \right) \tag{8}$$
So, we need to determine which is the greater leakiness comparing both fractions. But for reasonable parameters found in the literature[7], one might notice that \(\eta_2 \) is the minimum value achieved.
This determines that for the system fully off, circuit 2 minimizes the production of protein. The same steps can be followed to find leakiness for the half active system, meaning that only one of the inputs is present, and circuit 2 continues as the better one. This is crucial, to the therapy, because it means that the second system causes the least damage outside the targeted area.
Minimize the out-target protein production is essential, but now that the conceptual design was chosen, we need to find promoters that will maximize the production in the ON state and turn on precisely. Not only that, but we need promoters that are orthogonal to the E. coli Nissle 1917, so our system won't be affected by the bacteria's inner machinery.
As shown before, we can have a graph of the target protein utilizing the following substitutions. This will be a graph with two inputs in the x and y-axis (Lactate and Oxygen), and one output in the z-axis (Target protein). This is interesting because parameters can be checked and how they affect the slope of the 3D graph, utilizing the directional derivative in the direction of the vector (1,-1) because it will be activated in the ON region seen in Figure 9:
Figure 9: Hypothetical graph of GFP depending on the two inputs of the circuit.
If we were to take the directional derivative of GFP depending on the concentration of Lactate and \(O_2\), we would note that it depends on the Hill's coefficients. It is not actually necessary to do so, as we know the derivative multiplies the exponent. But you would need to take:
$$ D_{(1,-1)} = \boldsymbol{\nabla} GFP \boldsymbol{\cdot} (1,-1) \tag{9}$$
We conclude that the Hill's coefficients are very important for the slope, as expected because as n goes to infinity, we have a step function. (See figure 5).
Following the results, promoters need to be chosen in a way that they work with the second design, are orthogonal to the bacteria, have strong Y parameters and minimal y parameters and high Hill coefficients. The inputs promoters won't change, but the support promoters can be picked from Voigt et al. where he gives all constants fitted with a Hill function.
The deterministic mathematical model gave us important insight into the leakiness and helped us to go on with the design, also expliciting the importance of each parameter, being decisive to an optimal promoter choice and improving our circuit.
AND gate validation
Let's say we order the sequences, make the (very expensive) experimentations, just to see that our AND gate doesn't work, not an experimental mistake, but a logic one. Well, we believed it would function expected from its topology, but it isn't enough to guarantee a global behavior. How can we make our logic more robust, also helping to work out expected values? Of course, this is a job for the model team!
Let's begin where left off, from the circuit building. From the last section, we know the chosen design was from circuit one, and the promoter choices were made mainly from Voigt et al, as we ended up with the final design shown below:
Figure 10: A representation of our circuit is shown above, for more information, enter our design page: https://2019.igem.org/Team:Amazonas-Brazil/Design
As this is the same design as explained earlier, the equations are the same, only changing the placeholder promoters with the chosen ones. So, we conclude that the same steps explained in equations 3 and 4 can be utilized to find a general solution for GFP expression as \(t \rightarrow \infty \). Utilizing this method of substitutions, we find an expression of GFP depending on Lactate level, Oxygen and parameters. Meaning we can plot it on a 3 axis cartesian coordinate system, where z = GFP \( (Lactate , O_{2} ) \) . This is already an awesome result, as the AND gate can be exressed as a 2 input and 1 output function, that is ( GFP: \( {\rm I\!R} ^ {+} \times {\rm I\!R} ^ {+} \rightarrow {\rm I\!R} ^ {+} \) ). Doing the necessary procedures, we arrive at the equation:
$$GFP(Lact, O_2) = \frac{1}{\delta_{GFP}} \left( y_{G} + \frac { \bar{Y}_{G} (\delta_{T} R T) ^ {n_{T}} } { K_T(\delta_{T} R T) ^ {n_{T}} + V^{n_{T} } } \right) \tag{10}$$
Where. T represents constants related to TetR, \(\bar{Y} \) is the maximal - minimal production, and:
$$R = \delta_{L} ^ {n_{L}} K_{S} + \mathcal {H}_{L} ^ {n_{S}} \tag{11.a} $$
$$ T = \delta_{O_2} ^ {n_{O_2}} K_{P} + \mathcal {H}_{O_2} ^ {n_{P}} \tag{11.b}$$
$$ V = 2 y_{T}RT + T\bar{Y}_{T}\delta_{L} ^ {n_{S}} K_{S} + R \bar{Y}_{T}\delta_{O_2} ^ {n_{P}}K_{P} \tag{11.c} $$
And \(\mathcal {H} \) is the Hill function associated with each input, the S and P subscripts represent Srpr and PhiF, respectively.
Even though this seems pretty complex, don't be afraid, it is a great thing that this system can be solved analytically, as we don't need to rely on software to approximate our results, we get exactly what the model predicts. Also, think of all of this confusion as numbers, the only variables are Lactate and Oxygen, which will be our x and y coordinates and GFP our z coordinate, giving us a surface in the \( {\rm I\!R} ^ {3} \).
By itself, this is already an excellent result for any AND gate, which cannot be obtained through software. To understand why, we need to observe that a numerical analysis would give us GFP depending on time, for specific concentrations of Lactate and Oxygen, leading us to decide which concentrations are considered a "1" or a "0" for the inputs. This wouldn't be the best for the "gray area", where the system is not truly ON nor OFF. Our analytical result, however, plots the output as a function of the inputs, exactly what every AND gate wants!
The first important aspect that we arrive from the model only, which is great for the therapy, is that there are no output oscillations for fixed inputs, which is great expression wise, as it maximizes production when needed. Of course, you might be wondering: "Why is it a good result if it only works on the limit when time goes to infinity? ". Now, answering the question: the model predicts that this aspect occurs in the limit, but for the real-life, this is not always the truth. For example, when calculating how long does it take to completely charge a capacitor, you get that it takes infinite time, which is not true. But the model predicts that 99,999% is charged in a measurable time, meaning that in reality, 100% will be charged. This problem occurs as a consequence of considering something quantized as continuous, in the example, electrons are quantized. In our case, however, the copies of proteins are quantized. So, we can actually measure how long is long enough for our system, and we guarantee that it won't take too long, as the time is multiplied by the degradation rate (equation 4), which tends to be a very small quantity.
Okay, great! We know there is a function GFP(Lactate, Oxygen), but this does not mean we have an AND gate, as this function could behave in many unexpected ways, not really guaranteeing an AND gate. So, what does it take? Taking a look at Figure 5 again, we see how a toggle should work, if there are the inputs, there is the output, very simple. Now in two dimensions, would be the same thing to have a rectangle shape, as it would only be active as there are the two necessary inputs to have an output. For our continuous case, we expect something like what is observed in Figure 8.
To prove that our function (equation 10) has the expected behavior, we need to use calculus' tools. The first important thing to observe is that only looking at equation 11.a, we see that whenever Lactate grows(with constant oxygen), R grows, which means GFP grows and the opposite happens with oxygen. In other words, GFP is monotonically increasing in the x-axis and monotonically decreasing in the y-axis. And the critical points:
$$ \boldsymbol{\nabla} GFP = \vec{0 } \tag{12}$$
The Gradient can then be found using the chain rule:
$$ \frac{dGFP}{dLact} = \frac{dGFP}{dR}\frac{dR}{Lact} \tag{13.a}$$
$$ \frac{dGFP}{dO_2} = \frac{dGFP}{dT}\frac{dT}{O_2} \tag{13.b}$$
Then, we notice that the only two critical points for each input of the gradient are 0 and the limit to infinity. This means that the two critical points for GFP are (0,0) and \( (Lact \rightarrow \infty, O_2 \rightarrow \infty ) \). Not only that, but we notice that the function is limited. In the end, we have a limited monotonical function whenever one concentration is constant, two critical points on infinity and zero. This implies that every part curve with one constant input is sigmoidal (S/inverse S shape). To picture this, imagine that our function is a 3D object, then cut a very thin part parallel to x or y-axis, this cut has an S shape, like the one observed in Figure 5. When we implement this aspect in both axis, we get something we should expect for an AND gate with two inputs, basically a 3D sigmoidal function, a lot like a plateau.
Let's take a minute to appreciate this result. The steps shown above demonstrate elegantly how models predict behaviors, in this case, validates an expected one, adding robustness to the logic, making the work itself more scientific. To summarize, we modelled the system with deterministic functions, found a way to solve it analytically for time \(\rightarrow \infty \), founding a function of GFP concerning time and then used calculus tools to predict the behavior of the GFP, showing the expected continuous AND gate behavior.
Tumor growth model
Our project strives to become a therapy, but as it is known, tumors can be very different from one another, making it unlikely to exist an efficient generalized treatment. Another aspect to take into account is that we can't test how efficient the therapy is since experimentation won't be done. Of course, then it is a job for the modelling team to give us some idea of how it would work when implemented in real life. Not only that, but this might also be a great tool for someone interested in utilizing our therapy since, from indirect data, parameters can be measured and applied in the model, showing how the therapy would work in the specific tumor chosen.
The chosen approach was to build a populational model of ordinary differential equations. This, of course, is already a simplification, since a more accurate one would also take into account the position of the elements, turning it into a partial differential equation, which might be pretty untreatable depending on what you have. For our purpose, a populational model is good enough, since we only want to have an idea of the system might work.
$$\frac{d[C]}{dt} = \alpha C - \beta \frac{D}{K + D} \tag{14.a}$$
$$\frac{d[D]}{dt} = \omega B - \delta D \tag{14.b}$$
$$\frac{d[B]}{dt} = rB\left(1 - \frac{B}{\zeta} \right)\tag{14.c}$$
Where all concentration are in [g/cm³], being C the tumor cell density, D the drug density and B the bacterial density. The first thing to note is that we chose an exponential growth and a term of Drug action on tumor cells, for the tumor and a logistic for the Bacteria, where \(\zeta \) is the carrying capacity. And for the Drug, we chose a production rate that depends on Bacterial density and a degradation term, \(\delta \). Of course, depending on the drug chosen, other kinetics should be considered, but we are assuming that the presence of drugs is enough to kill tumor cells, not considering more complex mechanisms. Also, an important thing to consider is that the \(\omega \) parameter is related to the production of drugs when the system is completely ON, in a certain degree, this can be tunable for a specific cancer, making so it won't require much on the bacteria while maximizing tumor death rate[8].
As we know from the introduction section, these equations form a dynamical system. The usual method is to put the system in matrix form, find its eigenvalues and eigenvectors, enabling the use of the phase portrait. But, it is easy to see that our system is nonlinear, mainly looking at equation 14.c, meaning we can't put it in matrix form. What we can do, however, is linearizing the system in of its stable points, where all rates are equal to zero. This can be achieved through the Jacobian matrix which is for a 3x3 system:
$$ J = \begin{pmatrix}
\frac{\partial u_1}{\partial x_1} & \frac{\partial u_1}{\partial x_2} & \frac{\partial u_1}{\partial x_3} \\
\frac{\partial u_2}{\partial x_1} & \frac{\partial u_2}{\partial x_2} & \frac{\partial u_2}{\partial x_3} \\
\frac{\partial u_3}{\partial x_1} & \frac{\partial u_3}{\partial x_2} & \frac{\partial u_3}{\partial x_3}
\end{pmatrix} $$
For our proposed system, the stable point of interest is the origin, as it indicates that there are no more tumor cells. This leaves us with the following Jacobian:
$$ J_{(0,0,0)} = \begin{pmatrix}
\alpha & \frac{- \beta }{K} & 0 \\
0 & - \delta & \omega \\
\frac{r}{\zeta} & 0 & r
\end{pmatrix} \tag{15}$$
Now, we need to find it's eigenvalues and eigenvectors, solving the equation:
$$ det( {\bf J} - \lambda {\bf I} ) = 0 \tag{16} $$
This gives us the equation, where \(\lambda \) are the eigenvalues:
$$ \lambda ^{3} - (\alpha - \delta + r)\lambda ^{2} - (\alpha \delta - r \alpha + r \delta) \lambda + \left(\alpha \delta r + \frac{\alpha \omega r}{K \zeta}\right) = 0 \tag{17} $$
Even though we have the equation, solving a cubic is pretty convoluted. But, we know we want a negative eigenvalue, meaning that the origin is a local sink for the evaluated direction. So we define the variables, where b,c and d are the coefficients of the cubic equation 17:
$$ f \equiv \frac{3c -b^2}{3} \tag{18.a}$$
$$ g \equiv \frac{2b^3 - 9bc + 27d}{27} \tag{18.b}$$
$$ h \equiv \frac{g^2}{4} + \frac{f^3}{27} \tag{18.c}$$
For h<0, we have that all roots are real, and this is something we need since this linearization approach does not hold up for complex eigenvalues. If all roots are real, the solution goes as follows:
$$ j = \left(\frac{g^2 - 4h}{4} \right)^\frac{1}{6} \tag{19.a}$$
$$ k = arc cos \left(-\frac{g}{2j^3} \right) \tag{19.b} $$
$$ M = cos (\frac{k}{3}) \tag{19.c} $$
$$ N = \sqrt{3} sin \left(\frac{k}{3} \right) \tag{19.d} $$
$$ P = -\frac{b}{3} \tag{19.e} $$
Then, the eigenvalues are:
$$ \lambda_{1} = 2j cos \left(\frac{k}{3} \right) - \frac{b}{3} \tag{20.a}$$
$$ \lambda_{2} = L(M+N) + P \tag{20.b}$$
$$ \lambda_{3} = L(M-N) + P \tag{20.c}$$
This linearization with negative eigenvalues guarantees that the origin is locally stable in a certain direction, not globally. If we wanted to prove the origin as globally stable, we would need to find a Lyapunov potential function, which for our case is impossible, since every point in the direction (0,0, C) (no treatment case) diverges from the origin, as shown in Figure 11.
Figure 11: Qualitative behavior expected from the system with negative eigenvalues for linearization around the origin (red dot).
Now suppose someone wants to apply our treatment but isn't sure if it will be effective for a certain kind of cancer. Testing the therapy directly could be very expensive and hard to do, and without an idea of effectiveness, it is not the ideal case. So, this model offers an alternative way to indirectly test it, by measuring parameters such as tumor growth rate, the effectiveness of the chosen drug, bacteria growth rate and also a tunable parameter \(\omega \). And this model could also be adapted if other kinetics were to be considered.
Conclusion
We were able to successfully help the design team with our deterministic model, choosing the topology that minimized leaking, which is important for every new cancer treatment to reduce out-target damage. Afterward, our model helped to validate mathematically our implemented logic, confirming in a more robust way that GFP concentration in the steady-state, really formed an AND gate with 2 inputs, Lactate and \(O_2 \). And, at last, we built a simple Tumor growth model, which can help the implementation of our proposed therapy in other labs, with a tunable parameter and a condition for tumor decreasing locally. To wrap this up, despite working without experimental data, we were able to model and come to observable conclusions, moreover, we arrived at conclusions that are impossible to be determined when working only with numerical results given by software with measured parameters. An example is the output as a function of the inputs, which strongly enhances the power to predict behaviors, showing how important analytical analysis is, even though it is many times left in the dust. Thank you for reading and keep modelling!
References:
[1] Li, J., Georgescu, M., Hyde, P., Mahalov, A., & Moustaoui, M. (2015). Regional-scale transport of air pollutants: impacts of Southern California emissions on Phoenix ground-level ozone concentrations. Atmos. Chem. Phys., 15, 9345–9360.
[2] Elowitz, M., Leibler, S. (2000). A synthetic oscillatory network of transcriptional regulators. Nature volume 403, 335–338.
[3] Gardener, T., Cantor, C., Collins, J. (2000). Construction of a genetic toggle switch in Escherichia coli. Nature 403, 339–342.
[4] Strogatz, S. (1994). Nonlinear dynamics and Chaos.
[5] Jaoudé, W., Chaves, M., Gouzé, J. (2011). A Theoretical Exploration of Birhythmicity in the p53-Mdm2 Network.PLOS ONE 6(2): e17075.
[6] Goutelle, S. Maurin, M., Rougier, F., Barbaut, X., Bourguignon, L., Ducher, M., Maire, P. (2008). The Hill equation: a review of its capabilities in pharmacological modelling. Fundamental & Clinical Pharmacology, 22(6), 633-648.
[7] Stanton, B., Nielsen, A., Tamsir, A., Clancy, K., Peterson, T., Voigt, C. (2014). Genomic Mining of Prokaryotic Repressors for Orthogonal Logic Gates. Nat Chem Biol. 10(2), 99–105.
[8] Ratajczyk, E., Ledzewicz, U., Leszczynski, M., Friedman, A. (2017). The role of TNF-\(\alpha \) inhibitor in glioma virotherapy: A mathematical model. Math Bioscience Eng. 14(1), 305-319.
