In vivo expression of drugs
Binary oscillation gene circuit
In order to detect whether the drug can be expressed in the S.invicta as expected, we simulated and predicted the designed gene circuit. In the design and construction of gene circuits, we refer to and adopt the binary oscillation
model. The following is divided into four parts to describe it in detail.
Background
In nature, many organisms sense their surroundings by releasing small chemical signals and responding accordingly. For example, bacteria detect cell density by releasing and respond the self-induce signals that trigger the expression of
specific genes at the appropriate time. This process of cell-to-cell communication through chemical signals is known as Quorum Sensing (QS) . In recent years, researchers have discovered the existence of binary oscillations in the
process of QS, in which cell density repeatedly switches between two steady states. This oscillation phenomenon is important for understanding the bacteria behavior and regulation of organisms.
Purpose
The significance of this project is to gain an in-depth understanding of the mechanism of the binary oscillation phenomenon in the process of group sensing, and to explore the influence of the phenomenon on the survival and adaptation
of organisms to the environment. By establishing a mathematical model to simulate and analyze the binary oscillation phenomenon, we can reveal the intrinsic law of the phenomenon and provide theoretical support for subsequent experimental
design and observation. In addition, this study also helps to expand our understanding and application of ODEs in dynamic systems, providing new ideas and methods for research in related fields.
Model Content
Overview
In this model, we will address the following questions:
1. To develop a suitable model of the dynamical system ODEs to simulate the phenomenon of binary oscillations induced by populations in organisms. The model should consider key factors such as cell density, self-inducer concentration,
and their interactions.
2. The accuracy and applicability of the proposed model will be verified by experiments. We will adjust the model parameters and structure according to the experimental results so that it can better reflect the actual situation.
3. Analyze the kinetic behavior of the constructed model to explore the mechanism of the binary oscillation phenomenon and its influence on the survival and adaptation of organisms to the environment. We will reveal the intrinsic law of
the phenomenon by combining numerical simulation and theoretical analysis.
The Specific Content of The Model
The complete and specific formulas in the following contents can be found in the tex file.
1.For the simplified oscillation model
Fig.5 The simplified oscillation model
The ODE equations are:
The amount of A can be obtained as:
The sufficient conditions for further oscillation are as follows: 1.To ensure that the degradation rates of the key proteins A, B and C are similar, 2.and the degradation rate is appropriate to ensure that the group sensitivity threshold
can be activated.
Before we compare
and ,here we add assumption that A,B and C have similar degrade rate,i.e.dA≈dB≈dC.Thus
we can take
Now can be approximated to
and here we can derive that
Then subtract ,as a result
The oscillation period can also be directly calculated. It can be seen that the lower the degradation rate, the longer the oscillation period. Properly increasing the degradation rate can slow down the oscillation period and make the oscillation
more obvious, but too much may not reach the group sensing threshold.
2.For gene circuits
Fig.6 the binary oscillating system genetic circuit
The complete ODE equations are listed:
The numerical simulation results are
Fig.7 Binary oscillation curve result 1
Fig.8 Binary oscillation curve result 2
Summary
In nature, many organisms sense their surroundings by releasing small chemical signals and responding accordingly. The oscillation phenomenon is a big step in understanding the behavior and regulation of organisms.
This oscillation model aims to deeply understand the mechanism of the binary oscillation phenomenon in group sensing and explore the influence of some parameters on oscillation.
In addition, the theoretical analysis of the model also helps to expand our understanding and application of ODEs in dynamic systems, providing new ideas and methods for research in related fields.
Our contributions can be listed as follows:
1. We developed a suitable model of the dynamical system ODEs to simulate the phenomenon of binary oscillations induced by populations. The model considers critical factors such as cell density, self-inducer concentration, and interactions.
2. We verify the accuracy and applicability of the proposed model by experiments. The parameters of the model can be adjusted to adapt to different oscillations.
3. We gave some theoretical analysis to reveal the intrinsic law of binary oscillation. We proved the conditions for the oscillation so that, under some simple assumptions, we can ensure that the oscillation occurs and calculate the period
and amplitude of the oscillation.
Digital circuit gene circuit
For the previously designed binary oscillation gene circuit, although the expected results can meet our requirements, it has encountered difficulties in real experiments, and the binary oscillation model system is too complex and is likely
to be counterproductive. Therefore, we construct a new set of digital circuit genetic circuits to simplify the binary oscillation model.
Model Content
Overview
In this model, we will address the following questions:
1. To develop a suitable model of the dynamical system ODEs to simulate the phenomenon of binary oscillations induced by populations in organisms. The model should consider key factors such as cell density, self-inducer concentration,
and their interactions.
2. The accuracy and applicability of the proposed model will be verified by experiments. We will adjust the model parameters and structure according to the experimental results so that it can better reflect the actual situation.
3. Analyse the kinetic behaviour of the constructed model to explore the mechanism of the binary oscillation phenomenon and its influence on the survival and adaptation of organisms to the environment. We will reveal the intrinsic law
of the phenomenon by combining numerical simulation and theoretical analysis.
Parameter Description
Table 1 Parameters of chemical equation in gene circuit of digital circuit
| Symbol |
Definition |
|
|
The synthesis rate of nutrients
|
|
Synthesis rate of drug 1
|
|
Synthesis rate of LM17
|
|
The synthesis rate of LasI
|
|
The synthesis rate of LM17-LasI complex formed by the combination of LM17 and LasI
|
|
The synthesis rate of the rate at which the complex LM17-LasI initiates the expression of drug 2
|
|
The synthesis rate of the rate at which the complex LM17-LasI initiates the expression of TraR ( W )
|
|
The synthesis rate of the rate at which the complex LM17-LasI initiates esaI expression
|
|
Synthesis rate of TraR ( W ) -esaI complex formed by TraR ( W ) and esaI
|
| 0 |
The synthesis rate of the rate at which the complex TraR ( W ) -esaI initiates the expression of the cleavage product
|
|
Degradation rate of nutrients
|
|
Degradation rate of drug 1 |
|
Degradation rate of LM |
|
Degradation rate of LasI |
|
The degradation rate of complex LM17-LasI |
|
Degradation rate of drug 2 |
|
The degradation rate of TraR ( W ) |
|
Degradation rate of esaI |
|
The degradation rate of TraR ( W ) -esaI complex |
| 0 |
Degradation rate of pyrolysate |
The Specific Content of The Model
1.Expression of single bacteria in vivo
Considering the reaction and influence between substances:
① Biological safety device part:
The expression of nutrients is controlled by the anaerobic promoter, that is, E.coli can secrete nutrients for its own survival under anaerobic or anoxic conditions ( in the intestine ), and E.coli will die due to lack of
nutrients under aerobic conditions.
② The core part of the line :
First of all, this line was also initiated by the anaerobic promoter at the beginning, expressing
LM17,
LasI , drug 1 ( protease inhibitor ), involving quorum sensing, so there are the following reactions :
LM17 and LasI are combined into a complex, denoted as [ LM17-LasI ].
Then, when LM17-LasI reaches the threshold required for quorum sensing, it will make the next promoter express and promote the production of drugs 2 ( toxic protein ), TraR ( W ) , esal.
TraR ( W ) and esaI are combined into a complex, denoted as TraR ( W ) -esaI .
When TraR ( W ) -esaI reaches the threshold required for quorum sensing, it will make the next promoter expression, so that the expression of the cleavage gene.
The basic equations are listed according to the above reactions and effects :
When LM17 is less than the quorum sensing threshold of LM17 ( = 29.4 ), = 0;
When LM17 is greater than LM17 quorum sensing threshold ( = 29.4 ), = kff6;
Similarly, and .
When esaI is less than the quorum sensing threshold of esaI ( = 10 ), 0 = 0;
When esaI is greater than the quorum sensing threshold of esaI ( = 10 ), 0 = kff10.
Model Simplification
Assuming that the material degradation rate is negligible, all the above ODE can be simplified as :
According to the simplified ode equation finally determined above, we use simbiology in MATLAB to simulate and obtain the following results ( source file see : newmodel2 one.sbproj ):
Fig.9 Single E.coli expression line in vivo
The parameter settings : all kf is 0.7nM/s，kr is 0.07nM/s.
Operation results :
Fig.10 Single E.coli expression line in vivo
It can be concluded that drug 1 and drug 2 are delayed expression and positively correlated.
2.In vivo expression of multiple ( a group ) bacteria
Next, we consider the overall expression of a group of E.coli and E.coli lysis genes after expression. The change of OD600 during E.coli lysis measured in other team’s project was referred to：
Fig.11 The figure from the website https://2020.igem.org/Team:SZ-SHD/Proof_Of_Concept for details
Because the Y axis is OD600, we need to convert the number of bacteria. First of all, OD600 is proportional to the concentration of bacteria (See the website
https://zhuanlan.zhihu.com/p/648714725 for details.) , we can calculate the concentration directly according to the 8th power multiplied by the order of magnitude
10, and then set the volume to 10.Then the number of bacteria is the corresponding OD600 value multiplied by the 9th power of 10. After the conversion is completed, the red curve in the above diagram is fitted with python : ( code
see OD600.py file ).
Equation are obtained:
The curve image is as follows:
Fig.12 This is the curve of the number of bacteria over time fitted according to the above items.
The Y axis of the function is normalized, because our number wants to be expressed as a percentage, so the Y value of the first data point is changed from 1.47 to 1, and the subsequent Y value is reduced proportionally.
Continue to use python to derive the normalized function (code see OD600.py file ), and get the change rate curve of the number of bacteria :
y’=-0.1430299*exp(-0.209*x)
Fig.13 Derivative curve of the curve of bacterial number changing with time
Substitute into MATLAB to modify the original modeling:
Fig.14 Multiple E.coli expression lines in vivo and overall accumulation
The formula for the number of bacteria is as follows :
Operation results :
It can be seen that the number of bacteria decreases rapidly after lysing gene expression, and finally decreases to about 30 % and tends to be stable.
Fig.15 The curve of the number of E.coli lysed over time in the genetic route we designed
Fig.16 The cumulative changes of drug 1 and drug 2 in all bacteria
The peak expression of drug 1 and drug 2 was around 11 s. At this time, the amount of drug 1 is about 70nM, and the amount of drug 2 is about 40nM.
Summary
By using MATLAB and Python to simulate the digital circuit gene circuit, we can find that the results obtained are consistent with our expected experimental results. First of all, for a single cell analysis, we can draw the conclusion
that drug 1 and drug 2 are delayed and positively correlated. In fact, we must ensure that the protease inhibitor is produced one step ahead of the toxic protein, so as to ensure that the toxic protein in the S.invicta can be
accumulated without interference. Then for multiple bacteria were analyzed, we refer to the results of other teams in previous years, according to the mathematical calculation of the amount of bacteria equation. Substituting it into
the original model, we can get new results about the total expression of the population in vivo. According to the results, the number of bacteria decreases rapidly after lysing gene expression, and finally tends to be stable, rather
than falling to almost nothing. At the same time, we can also determine the peak expression levels of drug 1 and drug 2, so as to confirm whether our design meets the limits of the survival of S.invicta themselves, and whether
it has a negative impact on biosafety, the environment, etc.
Drug transmission and infection
Epidemic Model
In order to explore the spread and infection of endotoxin in ant colonies, and to predict the range and effect of our drugs after being expressed in vivo after being put into the S.invicta population, we established an infectious
disease model for the spread of drugs. This project uses the most basic SIR model.
After control
When we feed the medicine to the S.invicta, the S.invicta will behave like an epidemic, spread the medicine within the population, and finally be cured (to meet our requirements, the medicine can successfully work in the
S.invicta is considered to be cured. Our model is based on the SIR (Susceptibility, Contagion, Recovery) epidemic model. Here are some basic properties and assumptions of our model.
Assumptions
1.The population growth of S.invicta follows the curve in growth modeling. (Refer to Ant Colony Growth Modeling File )
2.The number of S.invicta population has been in a stable stage.
3.All S.invicta are susceptible individuals.
4.In our modeling, the healer is that the drug can successfully act in the body ( reaching a certain concentration ) ( here set the product at the beginning of the line to increase. )
5. S.invicta colony, ant nest as a closed area.
6. S (t), I (t), R (t) are continuous derivable functions of time t.
7.There is no incubation period of infection, and the onset of infection occurs within a short period of time after infection (consistent with the modeling of the expression of alternative routes in
S.invicta ).
8.Other factors that may affect the experiment can be ignored.
Parameter Description
Table 2 Parameters in the epidemic model
| Symbol |
Definition |
| S |
Susceptible people, groups that may eat bait
|
| I |
The patient, the group that eats the bait and plays a role
|
| R |
Convalescents, groups that did not function or died or recovered after eating the bait.
|
| N |
Total number of S.invicta
|
| S(t) |
The number of healthy S.invicta
|
| I(t) |
the number of infected S.invicta
|
| R(t) |
the number of removed S.invicta
|
| β |
Sympathetic rate
|
| γ |
Removal rate
|
| S(0) |
Initial value of health number
|
| I(0) |
Initial value of the number of infections
|
| R(0) |
Initial value of removal number
|
The Specific Content of The Model
According to the basic model of infectious diseases and the characteristics expressed by the
S.invicta after being fed with drugs, the ODE equation can be listed:
According to the setting of SIR, the value of the removal rate γ needs to be set smaller, set to 0.001, and the infection rate β is set to 0.3.Among them, (S ) in the second formula is the formula in the growth curve. Taking (S ) = 0,
the relationship and result diagram of the growth, infection and rehabilitation of the
S.invicta population were obtained by MATLAB programming. The code details are shown in the files SIRcode.txt and sir _ model.m.
Result
According to the infectious disease model equation listed above, the results obtained by programming are as follows ( the following time units are days ) :
Fig.17 Standard SIR model
The horizontal axis in the figure is time, and the unit is day; the vertical axis is the proportion of the corresponding number, which is the percentage. It can be seen that when the number of days is about 25 days, the infection rate
reaches its peak. At this time, almost all individuals in the S.invicta population eat the bait and our drugs are successfully expressed in their bodies. Therefore, it can be proved that the drugs we feed are effective enough
to spread to the entire population, so that there are toxins in each individual, and ultimately may achieve the purpose of feeding the toxins back to the queen and killing the queen.
According to the results of the S.invicta population growth model, we substituted 25 %, 50 %, 75 % and 100 % of the final growth peak (=10^6)in the ant nest into the infectious disease model ( S ), and obtained the following results
:
25％
Fig.18 Results of the epidemic model at 25 % of the final growth peak
50％of the final growth peak：
Fig.19 Results of the epidemic model at 50 % of the final growth peak
75 % of the final growth peak：
Fig.20 Results of the epidemic model at 75 % of the final growth peak
100 % of the final growth peak：
Fig.21 Results of the epidemic model at 100 % of the final growth peak
It can be found that the trend of all images is similar, and the infection rate is almost close to the total number of individuals in the whole population after a certain period when the
S.invicta begins to grow, that is, almost all individuals have eaten the drug and expressed.
Parameter Sensitivity Analysis
1. Parameter β sensitivity analysis
In this section, we carry out the sensitivity analysis of the parameter β of the SIR model. The infection rate β determines the rate at the number of an infected S.invicta spreads disease to other S.invicta per day. We analyzed
different β values and γ values to understand the impact of these two parameters on the spread of the epidemic.
Firstly, we define the range of parameter β, and generate 10 different β values between 0.1 and 1.0 with 0.1 as the step size. These different β values will be used to analyze the behavior of the model under different propagation rates.
At the same time, we also define another parameter gamma, which represents the recovery rate, which is set to 0.001. In order to preserve the results of the model under different parameters, we create a result matrix and pre-allocate
enough space to preserve the state of the model at 1001 time points for each β value. The final output is the state diagram of the model at different time points, which intuitively shows the evolution of the SIR model under different
β values. Each β value is represented by different line colors, which facilitates the comparison of model behavior under different parameters. The results are as follows :
Fig.22 Parameter β sensitivity analysis
By observing the results, the following conclusions can be drawn:
①The initial stage of the spread of infectious diseases: In all cases, the curve began to show an upward trend, which represents an increase in the number of infected S.invicta. This is because only a very small number of S.invicta are infected in the initial conditions, and most of the S.invicta are susceptible to ant colonies.
②Propagation peak: With the increase of β value, the curve rises faster, the peak appears earlier, and the number of infected S.invicta at the peak is lower. This means that higher β values lead to faster propagation, but lower
peaks.
③Curve decline: After the propagation peak, the curve began to decline. This is because the number of infected ants in the colony gradually decreased, while the number of recovered ants in the colony gradually increased.
④The smaller the β value is, the longer the propagation duration is. It is observed that when the β value is small, the rise time of the curve is longer and the propagation duration is longer. This means that reducing the β value will
slow down the rate of infection, making the epidemic spread longer and the peak higher.
In summary, β sensitivity analysis shows the important influence of β parameter on the transmission rate and scale of infectious diseases. Increasing the value of β will lead to an increase in the speed of infection, but the peak value
is lower and the duration of the transmission is shorter ; reducing the β value will slow down the propagation speed, but the peak value is higher and the propagation duration is longer.
2. Parameter γ sensitivity analysis
For analyzing the sensitivity of γ, we first fix the β parameter to 0.3 and then define the range and step size of γ. We generate 10 different γ values from 0.0001 to 0.01 with a step size of 0.001. Still pre-allocated enough space to
be able to save the state of the model at 1001 time points for each γ value. The result map is still drawn to visually show the evolution of the SIR model under different γ values. Each γ value is represented by different line colors,
which is convenient to compare the model behavior under different parameters. The results are as follows:
Fig.23 Parameter γ sensitivity analysis
By observing the results of the map, we can draw the following conclusion:
①Curve trend: All curves showed a trend of rising first and then falling, which is a typical SIR model behavior, reflecting the evolution of the number of infected people over time.
②γ value and peak value: With the increase of γ value, the peak value of the curve becomes larger. This is because the high γ value may also cause the duration of the infected person to be longer, so before the peak of the infection rises
to the highest point, the infected person is still increasing, thus making the peak higher, and the time when the curve begins to decline is delayed.
③γ value and decline time: With the increase of γ value, the time when the curve begins to decline is later. This is because the high γ value leads to a longer duration of the infected person, and the time when the recovery individual
begins to exist will be later.
④γ value and curve parallelism: It is observed that with the increase of γ value, the speed of curve decline remains unchanged, and the curve is completely parallel. This shows that the effect of the γ parameter on the recovery rate of
the infected person may be linear, so the decline rate of the curve is similar under different γ values.
All in all, the higher the value of γ does not make the number of infections decline faster as we suspect. From our results, the parameter γ is not sensitive. We can consider the subsequent combination of β and γ for more in-depth analysis
experiments, focusing on the impact of β on the results, and finding the key factors in the regulation of infectious diseases.
Summary
According to our results, it can be found that when the number of S.invicta is 25 % peak, 50 % peak, 75 % peak, and 100 % peak, the trend is similar, and it can infect almost all S.invicta populations within a certain period.
This proves that our drugs can effectively carry out almost all S.invicta individuals, and play a role in successfully expressing toxins, which lays a very important foundation for the accumulation and transmission of toxins
to the queen through back-feeding behavior. Through the sensitivity analysis of the two parameters β and γ, it can be concluded that the influence of β on the model is significant, so it is sensitive, while the change of γ has little
effect on the change and influence of the model, so we think it is not sensitive enough. After that, we can focus on the parameter β and carry out more in-depth experiments and analysis.
Drug delivery and enrichment
Enrichment model
In order to explore the cumulative efficiency of toxin transmission between young ants, workers and queens, we constructed a toxin enrichment model to predict the final toxin content accumulated in queens, so as to determine whether our
drugs can eventually kill queens successfully.
In order to simplify the model, the worker ants and worker ants are not considered, and the worker ants are fed to the young ants.
Ce represents the concentration of toxic protein in workers, represents the transmission efficiency of toxic protein from larvae to workers, Iye
represents the transmission rate of toxic protein from larvae to workers, Ne represents the number of workers, and Iye is related to the number of larvae Ny. Since the workers of
S.invicta cannot take unlimited food, but will receive a limit of maximum food intake, the correlation between Iye and Ny cannot be a linear relationship similar to Iye = kNy, but an increasing function with decreasing derivative
function. We use log curve approximation, where , , are the weight and bias parameters used to fit the curve;
Cq represents the concentration of toxic protein in queens,
represents the transfer efficiency of toxic protein from larvae to workers, Ieq represents the transfer rate of toxic protein from larvae to
workers, and Nq represents the number of queens. Ieq is related to the number of workers, and the correlation and modeling method are similar to the process of toxin transfer from larvae to workers, using log curve approximation.
The number of young ants Ny and the number of workers Ne are determined by introducing the ant colony growth model. In the ant colony model, the ordinary differential equation of the number of eggs laid by the ant colony changing with
time is y (t), and the ordinary differential equation of the number of workers changing with time is u (t).
According to the literature research, we know that in the ant nest, the worker ants account for 72 % ~ 95 %, and the soldier ants account for 2 ~ 3 %. Although the proportion of young ants varies greatly and irregularly, it can be inferred
from the proportion of other ant species that the proportion is about 2 % ~ 25 %. We take the median proportion, Ny and Ne can have the simplified representation of the total number of ants in the ant colony model x (t).
Because there is no suitable data fitting, we temporarily set the transfer efficiency
and
to 0.5, and the parameters of Iye and Ieq are designed as
,
= 1;
,
= 2 ;
,
= 0, and the actual parameters may be different from this parameter. However, according to the derivation of the model itself, we know that these parameters
only affect the specific concentration and speed of enrichment, and have no effect on the overall trend of enrichment. Therefore, the qualitative conclusions drawn from the model have reference significance.
According to the model, we use the scipy package of Python to solve the differential equations, and use the matplotlib package to draw the toxic protein concentration in the workers after 120 days as follows.
Fig.24 The concentration of toxic protein in workers after 60 days
The toxins in workers and queens are plotted together as follows;
Fig.25 The toxic protein concentration in workers and queens after 60 days
It can be seen from the diagram that the toxins in workers and queens are on the rise, and the rate of toxins rising in queens is much higher than that of workers, which can ensure that after the toxic protein is poisoned preferentially,
the purpose of targeted killing of queens can be achieved. Because the proportion of young ants in the total number of ants in nature fluctuates greatly and has no obvious regularity [1], in order to explore whether the toxic protein
can still be enriched in the rear of the ants in the natural environment, we set the proportion of young ants to a random number between 2 % and 25 %. At this time, the concentration of toxic protein in the worker ants is as follows.
Fig.26 Under the proportion of random young ants, the concentration of toxic protein in workers after 60 days
The toxins in workers and queens are plotted together as follows:
Fig.27 Under the proportion of random young ants, the concentration of toxic proteins in workers and queens after 60 days.
It can be seen from the figure that the toxins can eventually be enriched along the direction of larvae- > workers- > queens, thus proving that our toxic proteins are robust in complex environments.
Summary
In this part, based on the ant colony growth model of cubic function, we established a toxic protein enrichment model based on ordinary differential equations by simplifying the cross-feeding relationship between ant colonies. The
concentration of toxic proteins in each ant species is mainly affected by the transmission efficiency and rate of toxins and the number of ant species that transmit toxic proteins to it. At the same time, we considered the effect
of the upper limit of feeding on the enrichment of toxic proteins. We used Python 's scipy package to solve the model we constructed, and observed that the toxic proteins in workers and queens were on the rise. The concentration
of toxic proteins in workers increased rapidly in the first 10 days, and then increased slowly at a constant rate. The enrichment rate of toxic proteins in queens is thousands of times that of workers, which ensures that the toxic
proteins are preferentially poisoned after killing the ants to achieve the purpose of targeted killing the ants. At the same time, by introducing random numbers, we simulate the complex situation in which the proportion of some
ant species in the total number of ant colonies fluctuates greatly in nature. The toxic protein can finally be enriched along the direction of larvae, workers and queens. Therefore, it is proved that our toxic protein enrichment
scheme has robustness in complex environments.
