Team:SZU-China/Epidemic Model
Epidemic Model
Introduction
We developed a mathematical epidemic model to predict the population dynamics of cockroaches infected by our Metarhizium anisopliae since we can not do real experiment with our gene-modify product. We then performed numerical simulations on the model and sensitivity analysis on some key parameters. By doing so, we expected to estimate the efficiency with lethal time after using our ‘cockroach terminator’ indoors, and analyzed the impact of each relative factors that influence our product so as to modify our design.
Natural condition
In natural condition indoors, due to environmental resistance like food, water and space, the population of cockroaches is more likely to follow a S-shaped growth curve (sigmoid growth curve), which can be formalized mathematically by logistic function $\dfrac {dN}{dt}=rN\left( 1-\dfrac {N}{K}\right)$ .After we research papers and government researches, we found that
After simulating the numerical solution (Figure 1), we can see that the cockroach population growing rapidly and reach to stationary state of 40 less than fifteen days. It’s necessary for as to control this worrisome infection.
Figure 1: The population growth curve of cockroaches in house from one to large amount.
After control
After we using M.anisopliae emulsifiable powder to control the cockroach population indoors, M.anisopliae will behave like an epidemic-like disease and spread out inside the population and finally kill the cockroaches. Our SID epidemic model with ordinary differential equation (ODE) was constructed based on SIR (Susceptible, Infectious, Recovered) epidemic model, following are some basic properties and assumptions of our model.
Assumptions
- 1. The population growth of cockroach follows sigmoid growth curve.
- 2. The population of cockroach has reached close to the stationary stage.
- 3. Ignore natural death rate in our system since stubborn vitality of cockroach.
- 4. Infectious individuals can not recover, and ignore they give birth.
- 5. Naturally all cockroaches are susceptible individuals, they can infect by M.anisopliae becoming infectious individuals.
- 6. The number of individual being infected in a contact between a susceptible and an infectious subject is determined by standard incidence rate $\dfrac {\beta IS}{S+I}$ .
- 7. The transition rate between infectious and dead is γ, its reciprocal ($\dfrac {1}{\gamma}$) determines the average infectious period, which is estimate by experiment data.
- 8. Other factors that may affect the experiment are negligible.
Figure 2: Diagram of epidemic model. Showing the relation between each group of individual.
Parameters
Table 1: Parameters and variables for the model.
| Parameter | Value | Meaning | Source |
|---|---|---|---|
| β | 0.775 | transmission rate, which is the probability of getting the infection in a contact between susceptible and an infectious | [2] [3] |
| γ | 1/7 | mortality, which is the the transition rate between I and D, its reciprocal (1/γ) determines the average infectious period | B |
| r | 0.61 d-1 | growth rate | [1] |
| K | 40 | carrying capacity | A |
| S(t) | variable | the number of susceptible individuals over time | - |
| I(t) | variable | the number of infectious individuals over time | - |
| D(t) | variable | the number of dead individuals over time | - |
| N(t) | S(t)+I(t)+D(t) | population size | - |
| S(0) | 35 | the initial number of susceptible individuals | A |
| I(0) | 3 | the initial number of infectious individuals | - |
| D(0) | 0 | the initial number of dead individuals | - |
|
[]See in the reference code at the bottom AGovernment researches from Shenzhen Center for Disease Control and prevention BEstimate by experimental data |
|||
By the above relations between each group of number, the system based on assumption described above can be expressed by the following set of ordinary differential equations: \begin{align*} &\dfrac {dS}{dt}=rN\left( 1-\dfrac {N}{K}\right) -\dfrac {\beta IS}{S+I} \\ \\ &\dfrac {dI}{dt}=\dfrac {\beta IS}{S+I}-\gamma I \\ \\ &\dfrac {dD}{dt}=\gamma I \\ \\ &\; N=S+I+D \ \end{align*} This system is non-linear, and the analytic solution does not exist, but we can compute the numerical solution. We used improved Euler method to estimate the discrete solution in one step size to simulate the change of number with time in days (see the results). And Runge-Kutta methods solution for variable-step by MATLAB to constructed the sensitivity analysis.
Results
The following curves (Figure 3) show dynamics of number change of each kinds of individuals.
We see that the infectious individuals grow fast before first 7 day, and then began to drop. The total number of cockroaches continuously going down.
We specify the median lethal time (LT50), which means 50% of cockroaches dead, where in this condition is about 8 day, and lethal time 80% (LT80) is about 15 day. This result is reasonable with our expectation, it proved that our product indeed control the total number of cockroaches efficiently in doors.
Figure 3: The epidemic model simulating change of each group of individual within 40 days.
Sensitivity Analysis
We use sensitivity analysis to analyze the impacts of some important parameter values (β, γ) on our model outcomes lethal time 50% and 80% (LT50 and LT80), in order to figure out what we should need to improve in future, and which parameters we need to be more certain.
We ran our model with transmission rate (β) and mortality (γ) ranging from 0 to 1 and plotted these values against the time in LT50 and LT80. What’s more, we introduced sensitivity coefficient (△t*p/t *△p ratio of the relative change of the lethal time to the relative change of the parameter) into consideration, increasing our parameters 10 percent, and find the relative change of the lethal time. The figures below show the changing tendency of lethal time with respect to each parameter and the sensitivity coefficient in ten percent up.
Figure 4: Change of the lethal time when transmission rate (β) ranging from 0 to 1. Blake dot represent our parameter used in Table 1.
Figure 5: Change of the lethal time when mortality (γ) ranging from 0 to 1. Upper x-axis is 1/γ, indicating the average mortality day. Blake dot represent our parameter used in Table 1, and dash line showing when lager than 0.8 because ODE system can not find solution in our interval.
We can conclude the intuitive results after scanning either parameter from 0 to 1. For transmission rate (Figure 4), as beta increasing, both LTD50 and LT80 going down. However, the change of lethal time become more smoothly when beta closer to one. The black dot in curves indicates our simulating parameter value 0.775 in epidemic model, which can consider as a good approximation and high efficiency that we do not need to increase so much.
As for mortality (Figure 5), lethal time will increase dramatically when gamma gets closer to zero or one (when gamma closer to one, the ODE system can not find real solution in our interval), but change smoothly in a low region in the middle. That’s sound great for as since the black dot, indicating our simulating parameter value 1/7 in epidemic model, lines near the middle. So, if we can control to increase mortality a little bite, it can make a difference. This is just what we do in our project. The functional genes we used may decrease lethal time into nearly stationary state.
Figure 6: Sensitivity coefficient of transmission rate (β) and mortality (γ) to LT50 and LT80.
Moreover, the sensitivity coefficient in Figure 6 shows the coefficient is less than -0.8, which indicate a low sensitivity of beta and gamma, while gamma has higher sensitivity compared with beta. In conclusion, the parameters we estimated previously did not greatly affect the lethal time of our model, we can say our parameters are well approximate.
Discussion
By our improved epidemic model, we successfully predict the population dynamics of cockroaches infected by our Metarhizium anisopliae indoors. Estimate the key parameters beta to be 0.775, gamma to be 1/7, can be considered as a good approximation not only we estimate them from research papers and experimental data, but also we constructed sensitivity analysis to these parameters and showed a low sensitivity. These results give us a clear relation about each factor and to direct our project design that gamma is the most important and easy factor we need to improved and make a difference. These results give us a clear relation about each factor and to direct our project design. Numerical simulations on the model showed our product kill cockroaches indoors with a lethal time 50% (LT50) about 8 day, and lethal time 80% (LT80) about 15 day. Thanks to the assumption, we can simplify our model, however in future, we will consider more factor that may influence our prediction to make our model more real.
References
[1] D. M. Müller-Graf, Jobet, E., Cloarec, A., Rivault, C., Baalen, M. V., & Morand, S. (2001). Population dynamics of host-parasite interactions in a cockroach-oxyuroid system. Oikos, 95(3), 431–440.
[2] Campaña, Ana & Funderburk, Karen & Kaur, Amandeep & Puente, Patricia & Ríos-Soto, Karen. (2017). A Household Model of Cockroach Infestation and Its Effects on Atopic Asthma Symptoms.
[3] Ni, Yu-Ting, et al. (2014). Metarhizium anisopliae as a Potential Microbial Agent for Managing the Brown-banded Cockroach, Supella longipalpa (F.).
[4] Hu, Z., Liu, S., & Wang, H. (2008). Backward bifurcation of an epidemic model with standard incidence rate and treatment rate. Nonlinear Analysis Real World Applications, 9(5), 2302-2312.