Cytomegalovirus dynamics model with random behavior

1.   Introduction

Cytomegalovirus (CMV) is a kind of herpesvirus DNA virus. The infection rate of adults in China is greater than 95%. CMV infection is mostly acquired in childhood. CMV is usually a recessive infection with no obvious symptoms in adults. The invasion of CMV into organs and systems can cause diseases, especially in infants, which can lead to body defects or death in severe cases. CMV can be transmitted through intrauterine infection of fetuses, breastfeeding of infected infants, blood transfusion or organ transplantation ^[1,2,3]. The treatment and prevention of CMV infection had always been an important research topic. At present, medical workers focus on the study of intensive treatment for infants infected with CMV and the prevention of vertical infection between mother and infant. According to relevant regulations, the infected CMV breast milk cannot be fed directly. It needs to be fed after the CMV is killed by breast milk handling. However, in real life, some CMV infected mothers do not handle breast milk and feed directly, which does not lead to infant infection. Whether or not to breast milk infected with CMV needs to be completely processed before fed has been a controversial issue ^[4,5]. The clinical data of CMV infection and breast milk treatment in children were analyzed. The analysis results show that there is a negative correlation between CMV infection and breast milk handling. This further confirms the important value of breast milk handling in preventing infant infection ^[6,7]. For children infected with CMV, drug treatment is one of the preferred methods, and the choice of medication regimen will affect the treatment to a certain extent. Through the experimental data analyzed of different medication regimens for CMV children, it is found that intravenous and oral combined medication have a better therapeutic effect ^[8]. Because of the interference of uncertainties such as treatment and breast milk handling, it is impossible to determine the information of medication and feeding through a large number of experiments and it is difficult to understand many problems of treatment and prevention of children infected with CMV. Using the viral dynamics model, the deviation between dosage and breast milk handling in predicting CMV infection in children is discussed theoretically, which solves the problem of inconvenient observation in clinical experiments and provides theoretical support for formulating scientific defensive measures. At present, there are a few studies on the dynamic model of human cytomegalovirus. However, the results of the viral dynamics model and antiviral treatment model are numerous, such as the model of hepatitis B, AIDS and other viral infectious diseases. These results provide a theoretical basis for the study of cytomegalovirus dynamics model ^[9,10]. Because of the uncertainty of people's living habits and medical compliance, it is practical to use the stochastic viral dynamics model to study the changing trend of viruses and control strategies. Mckendrick et al. found that there were uncertainties in human contact rate, cure rate and mortality rate in the study of infectious diseases, and the random method was used to describe the transmission process of infectious diseases, and an infectious disease model was established. The influence of the disturbance of various factors on the transmission law of diseases was studied ^[11]. Based on the data simulation and mechanism analysis of the number of SARS infections and patients, a random SARS epidemic model was constructed. It was found that the change of infection rate was the most important factor affecting the spread of SARS ^[12,13]. Under certain conditions, the model can predict the epidemic situation and random fluctuation trend as well. Zhang Lifeng ^[14] compared the difference between deterministic and stochastic models of infectious diseases and found that the two models were based on the basic reproductive number to judge the epidemic of infectious diseases. The infection rate and cure rate can be simulated accurately by counting the numbers of infected and susceptible persons. The estimation method of the basic regeneration number was given in ^[14] and ^[15]. Xia Peiyan's group studied the infection of HTLV-1 and HIV-1 virus. They considered the replication and death of the virus were interfered by factors such as environment and virus itself and considered that the virus was interfered by white noise. They established a virus dynamics model and studied the existence, ergodicity and asymptotic of the solution of the model. The results showed that the virus will always exist in the organism and fluctuate near the equilibrium point of the deterministic model. The magnitude of the fluctuation depends on the intensity of the disturbance ^[16,17,18].

By considering the mechanism of CMV infection and the existing virus models, the main contributions are as follows: 1) investigate the random phenomena in treatment and breast milk handling, 2) establish a stochastic cytomegalovirus dynamics model, and 3) investigate the existence, uniqueness and ergodicity of the model solution. The effects of treatment and feeding on children with infection are analyzed by numerical simulation.

In the second part of the paper, a cytomegalovirus dynamic model with random behavior is established according to the transmission mechanism of the virus and the mode of CMV infection. In the third part, the global existence and uniqueness, asymptotic behavior, ergodicity and extinction of the solution of the model are studied by using the stochastic differential equation theory. The fourth part uses the numerical simulation method to analyze the effect of parameter perturbation on the variation law of model solution and to verify the feasibility of the results of the theorem. In the fifth part, suggestions for the prevention of CMV infection are provided according to the results of the study.

2.   Model establishment

CMV is a DNA virus that can only grow in living cells. When the cell is infected, the size of the nuclei increases. CMV replicates in some undetermined monocytes such as mononuclear phagocytes and immune cells, we record the number of these healthy cells as $T$. CMV can cause damage to healthy cells, and we record the number of such infected cells as $I$. The number of virus that CMV will replicate is recorded $V$. Long-term infection can damage tissue, animal experiments ^[19] have verified the process of CMV infection in other cells. There is no specific treatment for CMV infection so far, so only symptomatic treatment can be taken. For children with severe symptoms of CMV, antiviral drugs can reduce CMV infection to healthy cells and reduce the number of infected cells effectively. As for the breastfeeding, unknown mothers infected with CMV without clinical symptoms and known mothers infected with CMV without breastfeeding directly can also increase the number of CMV in a baby and lead to a certain risk of infection ^[20]. According to the developmental characteristics of the newborn, the role of the infant's immune system against CMV infection is not considered. Only breastfeeding and drug therapy are considered. Assuming that breast milk contains $\gamma$ times as many CMVs as infants, which feds into infants directly, thus reducing the mortality rate of CMV in infants, which is recorded as ${b_1} = b - \gamma$. Based on the mechanism of replication of Nowak virus dynamics model ^[21,22], the deterministic cytomegalovirus dynamics model is established as follows

where $T$, $I$ and $V$ respectively represent the number of healthy cells, giant cells and free viruses at time $t$. $A$ is the number of growing healthy cells; $\beta$ is the infection rate of viruses infecting healthy cells; $\rho$ is the rate of viruses releasing from giant cells; $b$ is the natural mortality rate of healthy cells and infected cells; $\alpha$ is the reduced rate of infected cells under treatment.

The basic regeneration number of Model (1) can be obtained. We get

The virus-free equilibrium point is $E_{0} = (A / b, 0, 0)$ and the virus equilibrium point is ${E^*} = ({T^*}, {I^*}, {V^*}) = (\frac{{{b_1}(b + \alpha)}}{{\beta \rho }}, \frac{{\beta \rho A - (b + \alpha)b{b_1}}}{{\beta \rho (b + \alpha)}}, \frac{{A\rho }}{{(b + \alpha){b_1}}}(1 - \frac{1}{{{R_0}}}))$.

From the method in ^[23], we can get that the virus-free equilibrium point of Model (1) was globally asymptotically stable when ${R_0} < 1$, and the virus-free equilibrium point of Model (1) was globally asymptotically stable when ${R_0} > 1$.

Many viruses that invade the body produce antibodies and immune lymphocytes that limit virus replication, but they are not resistant to the activation of endogenous latent viruses and exogenous infection by different strains of viruses. This leads to the uncertainty of recurrent infection or recurrence after cure, resulting in uncontrollable virus phenomenon, known as "viral blip" phenomenon. For example, cytomegalovirus (CMV), human immunodeficiency virus (HIV), hepatitis B virus (HBV) and so on. There are many uncertainties in the infection and treatment of such diseases, which lead to this phenomenon ^[24,25,26]. In ^[27], Zhang Wenjing et al. studied a certain type of HIV anti-oxidation treatment model and immune model, and also found the existence of this phenomenon, and provided the conditions for the existence of the virus spot. Numerical simulations have shown that saturated infection patterns contribute to recurrent infections. Drug is one of the common means to control virus, but treatment strategy has certain influence on virus suppression. In ^[28], Wang Shaoli et al. studied the mathematical model of a class of single-strain and multi-strain virus infection, and found that the effect of different treatments may lead to the competitive rejection of virus strains, and proposed the suggestion of combined treatment.

For those infected with CMV immunodeficiency, CMV infection may be more serious, so that the body's immune suppression is unable to resist the role of the virus. Because the infant's immune system is not robust, it cannot use autoimmunity against CMV. After CMV infection, viral replication is sufficient to cause damage to target organs or tissues requiring treatment. If no target organ damage is caused, the virus is latent or incomplete. Effective measures should be taken if the child receives an exogenous CMV infusion into the foot to infect damaged organs and tissues. For example, breast milk was fed with or without infected CMV. Due to the randomness of CMV treatment in breast milk in real life, the amount of virus imported into infants is uncertain, which may cause infant CMV infection or agammaavate infection. For infected infants, uncertainties such as the choice of medication mode, drug leakage and insufficient dosage of medication can also occur in the treatment, which cannot achieve the desired therapeutic effect. These random phenomena can influence with the number of healthy cells $T$, the number of infected cells $I$ and the mortality rate $b$ of CMV in Model (1). The uncertainties will cause disturbance of Model (1) parameters. Assuming that

The following stochastic cytomegalovirus dynamics model is obtained from Model (1)

Let be a complete probability space and is the$\sigma$-algebra. In this probability space, a $\sigma$-algebraic stream ${\left\{ {{F_t}} \right\}_{t \geqslant 0}}$ is defined, which satisfies the usual conditions.

(1) For all$0 \leqslant s < t < \infty$,;

(2) Right continuity: for all$t \geqslant 0$, .

Moreover, the ${B_i}(t), i = 1, 2, 3$ in Model (2) belongs to the independent Brownian movement in this probability space, and ${\delta _i}$ is the wave force of the Brownian movement. In terms of biological significance, Model (2) should be investigated under the condition $R_ + ^3 = \left\{ {(T, I, V)} \right.$$\left| {\left. {T > 0, I > 0, V > 0} \right\}} \right.$ and $(\Omega, F, P)$.

3.   Main research results

3.1.   Existence and uniqueness of the global positive solution

Theorem 1. For any initial value $(T(0), I(0), V(0)) \in R_ + ^3$, there is a unique positive solution $(T(t), I(t), V(t))$ of System (2) on $t \geqslant 0$ and the solution will remain in $R_ + ^3$with probability one, namely, $(T(t), I(t), V(t)) \in R_ + ^3$ for all $t \geqslant 0$ almost surely (a.s.)

Proof. Since the coefficients of System (2) are locally Lipschitz continuous, then for any initial value $(T(0), I(0), V(0)) \in R_ + ^3$ there is a unique local solution $(T(t), I(t), V(t))$ on $t \in [0, {\tau _e})$, where ${\tau _e}$ denotes the explosion time ^[29].

To proof that the solution exists globally, we only need to verify that ${\tau _e} = \infty$ a.s.

Let ${m_0} \geqslant 0$ be sufficiently large and the values of $T(0), I(0), V(0)$ all lie within the interval $[m_0^{-1}, {m_0}]$. For each integer $m \geqslant {m_0}$, the stopping time is

where $\inf \varphi = \infty$ ($\varphi$ denotes the empty set). It can be seen that ${\tau _m}$ increases as $m \to \infty$. Let ${\tau _\infty } = \mathop {\lim }\limits_{m \to \infty } {\tau _m}$, whence $\tau {}_\infty \leqslant {\tau _e}$ a.s. If it can be proved that ${\tau _\infty } = \infty$ a.s. is true, then ${\tau _e} = \infty$ a.s. and $(T(t), I(t), V(t)$$(t \geqslant 0)) \in R_ + ^3$ a.s. If this assertion is false, then there exists a pair of constants $\bar t > 0$ and such that

Thus there is an integer ${m_1} \geqslant {m_0}$, so that $m \geqslant {m_1}$

Define a ${C^2}$-function $V$: $R_ + ^3 \to [0, + \infty)$

where ${c_1}, {c_2}$ are positive constants to be determined later.

Function $V(T, I, V)$ is nonnegative. Applying $It\hat o's$ formula to $V$, we get

where

Choose

such that

Then

Using Formula (6) to calculate Formula (4) obtainable

Integrating (7) from 0 to $t$ at both ends and then take the expectation

Then

Let ${\Omega _m} = \left\{ {{\tau _m}} \right. \leqslant \left. {\bar t} \right\}$. According to the definition of stopping time and Formula (3), for each $\omega \in {\Omega _m}$, there is $T({\tau _m}, \omega), I({\tau _m}, \omega)$ and $V({\tau _m}, \omega)$ at least one of them equal to $m$ or $1/m$, so

So we can see from Formula (8)

where ${J_{{\Omega _m}}}$ is the indicator function of ${\Omega _m}$.

Let $m \to \infty$, then$V(T(0), I(0), V(0)) + N\bar t > \infty$ is a contradiction. Thus${\tau _\infty } = \infty$. This means that $(T(t), I(t), V(t))$ with probability 1 does not produce blasting in a limited time period. The proof is complete.

3.2.   Asymptotic behavior of solutions

Theorem 2. If ${R_0} > 1$, and $\delta _1^2 < \frac{3}{2}b$, $\delta _2^2 < 2(b + \alpha)$, $\delta _3^2 < 2{b_1}$ are valid, then the solution of Model (2) has the following properties.

where

${E^*} = ({T^*}, {I^*}, {V^*})$ is the equilibrium point of Model (1).

Proof. Define a ${C^2}$function $V:R_ + ^3 \to {R_ + }$ as follows

Applying $It\hat o's$ formula to $V$, we get

where

When $a = \frac{{{{(2b + \alpha)}^2}}}{{b + \alpha - \frac{1}{2}\delta _2^2}}$, Formula (10) can be changed to

From Formula (11), and substituting $a$ into Formula (9), we can get

Integrating (12) from 0 to $t$ at both ends and then take the expectation

where

Formula (13) deformation is available

$E\int_0^t {\left\{ (\right.}$$\frac{3}{4}b - \frac{1}{2}\delta _1^2){(T - \frac{{3b}}{{3b - 2\delta _1^2}}{T^*})^2}$

Divide the two ends of Eq (14) by $t$ and let $t \to \infty$, It can be obtained

where $M = \min \left\{ {\frac{3}{4}b - \frac{1}{2}\delta _1^2, b + \alpha - \frac{1}{2}\delta _2^2, {b_1} - \frac{1}{2}\delta _3^2} \right\}$ is not negative, so the theorem holds. The proof is complete.

The results of Theorem 2 show the existence of perturbation and change the law of solution of the original Model (1). Although the stochastic perturbation makes Model (2) have no definite positive equilibrium point, its solution oscillates around a fixed point ${P^*} = (\frac{{3b}}{{3b - 2\delta _1^2}}{T^*}, \frac{{2{{(b + \alpha)}^2} + \delta _2^2({{(2b + \alpha)}^2} - (b + \alpha))}}{{b + \alpha - \frac{1}{2}\delta _2^2}}{I^*}$ at infinity in the sense of mean value. The amplitude of vibration is no larger than $K / M$ and it is related to the size of $\delta _i^2$. When there is no disturbance, which is $\delta _i^2 = 0$, there is a positive equilibrium point in the Model (2), and then ${P^*} = {E^*}$returns to the original equilibrium point.

3.3.   Ergodicity of solutions

Let $X(t)$ be a homogeneous Markov process in ${E^l}$(${E^l}$ represents $l$-dimensional Euclidean space), and it can be described by the following stochastic differential equation

Moreover, the ${B_i}(t), i = 1, 2, 3$ in Model (2) belongs to the independent Brownian movement in this probability space. The diffusion matrix is defined as follows

Lemma 1. (ergodicity theorem) ^[30] If there exists a bounded domain $U \subset {E^l}$ with regular boundary $\Gamma$ then:

${A_1}$ : In domain A and some neighborhoods, the minimum eigenvalue of the diffusion matrix is nonzero.

${A_2}$ : When $x \in {E_l}\backslash U$, the average time $\tau$ of the trajectory path from point $x$ to set $U$ is limited, and there is $\mathop {\sup }\limits_{x \in Q} {E_x}\tau < \infty$ for each compact subset $Q \subset {E_l}$.

Then the Markov process $X(t)$ of Eq (15) has a stationary (unchanged) distribution $\mu (\cdot)$.

Theorem 3. Assuming that $R_0^S > 1$, then Model (2) has a unique stationary distribution $\mu (\cdot)$ and it has the ergodic property.

where

Proof. Theorem 1 shows that the solution of Model (2) is existence and uniqueness.

The diffusion matrix of Model (2) is given by

Choose $W = \mathop {\min }\limits_{(T, I, V) \in {{\bar D}_\delta } \subset R_ + ^3} \left\{ {\delta _1^2{T^2}, \delta _2^2{I^2}, \delta _3^2{V^2}} \right\}$, , ${\bar D_\delta }$is a bounded subset of $R_ + ^3$, we have

Then the condition ${A_1}$ in Lemma1 holds.

Define

where ${c_1}$ and ${c_2}$ are positive constants to be determined later.

Applying $It\hat o's$ formula to $V$, we get

Let

Then

Formula (16) can be simplified to

where

Because $R_0^S > 1$, there is $\lambda > 0$.

For a sufficiently small constant $h$, take a positive number H and define a ${C^2}$ function $\bar V$:$R_ + ^3 \to {R_ + }$ as follows

Easy to verify

where $\varepsilon > 0$ is a sufficiently small constant

In addition, $\bar V(T, I, V)$ is a continuous function. Finding partial derivatives for $\bar V(T, I, V)$

Rearranging Eq (19)

Substitute Formula (20) into the first equation formula of (19), we have

There is a positive real root ${T_0}$, which is substituted by Formula (20) to obtain ${I_0}$ and ${V_0}$. So we know that $\bar V(T, I, V)$ has a minimum point $({T_0}, {I_0}, {V_0})$ in $R_ + ^3$.

Modified (18), define a nonnegative ${C^2}$ function $\hat V:R_ + ^3 \to {R_ + }$

Applying $It\hat o's$ formula to ${V_2}$、${V_3}$、${V_4}$, we get

where

Formulas (17), Formulas (21), Formulas (22) and Formulas (23) are available. We get

In order to prove that the condition ${A_2}$ in Lemma 1 holds on $D_\varepsilon ^C$ = $R_ + ^3\backslash {D_\varepsilon }$, we rewrite set $D_\varepsilon ^C$. The first six regions are as follows.

It can be found that $D_\varepsilon ^C = {D_1} \cup \cdots \cup {D_6}$, and the following formula can be verified in ${D_i}(i = 1, 2, 3, 4, 5, 6)$.

Case 1. If, for sufficiently small, by Formula (24), we have

where

is a finite number.

For all $(T, I, V) \in {D_1}$, by Formula (26) knows Formula (25) is established.

Case 2. If, for sufficiently small, by Formula (24), we have

where $C$ has an expression in the form of Formula (27), which is still A finite number in${D_2}$.

For all $(T, I, V) \in {D_2}$, according to Formula (28), we have deduced that Formula (25) is established.

Case 3. If, Choose positive by Formula (17), and for we have. For sufficiently small, by Formula (24), we have

where

is a finite number.

For all $(T, I, V) \in {D_3}$, according to Formula (29), we have deduced that Formula (25) is established.

Case 4. If, for sufficiently small, by Formula (24), we have

where

is a finite number.

For all $(T, I, V) \in {D_4}$, according to Formula (30), we have deduced that Formula (25) is established.

Case 5. If, for sufficiently small, by Formula (24), we have

where

is a finite number.

For all $(T, I, V) \in {D_5}$, according to Formula (31), we have deduced that Formula (25) is established.

Case 6. If, for sufficiently small, by Formula (24), we have

where

is a finite number.

For all $(T, I, V) \in {D_6}$, according to Formula (32), we have deduced that Formula (25) is established.

Summing up the evidence, for all $(T, I, V) \in R_ + ^3\backslash {D_\varepsilon }$, we have deduced that Formula (25) is established. Thus Model (2) is ergodic and has invariant distribution. The proof is complete.

Theorem 3 shows that the virus will eventually remain in a constant state for a sufficient time period.

3.4.   Extinction of solutions

If $f$ is an integrable function on $[0, \infty)$, define $\left\langle {f(t)} \right\rangle = \frac{{\int_0^t {f(s){\rm{d}}s} }}{{\rm{t}}}$.

Theorem 4. If $\hat R_0^s = \frac{{\rho \beta A}}{{b(b + \alpha)({b_1} + \frac{{\delta _2^2\delta _3^2}}{{2(\delta _2^2 + \delta _3^2)}})}} < 1$ holds, then the solution of Model (2) has $\mathop {\lim }\limits_{t \to \infty } \left\langle {T(t)} \right\rangle = \frac{A}{b}$, $\mathop {\lim }\limits_{t \to \infty } I(t) = 0$, $\mathop {\lim }\limits_{t \to \infty } V(t) = 0$ a.s.

Proof. Define the Lyapunov function

Applying $It\hat o's$ formula to ${U_2}$, we can get

where

From the first equation of Model (2), we can get

Integrating it from 0 to $t$, by Theorem 1 of ^[31], we can get

By Formulas (35) and $\hat R_0^s = \frac{{\rho \beta A}}{{b(b + \alpha)({b_1} + \frac{{\delta _2^2\delta _3^2}}{{2(\delta _2^2 + \delta _3^2)}})}} < 1$, Formulas (34) scaling has

Substitute Formula (36) into Formula (34), integrating Formula (34) from 0 to $t$, by Theorem 1 of ^[32], we can get

Thus we can see that

Add up the three equations of Model (2), we have

Integral to it, divide both sides by $t$ at the same time, we can get

From Formula (37) and Formula (38), we can have

The proof is complete.

4.   Numerical simulation

In order to understand the effect of the uncertainty of treatment and breast milk handling in killing CMV on the severity of CMV infection. Taking $\alpha$ and ${b_1}$ as the research objects, A and B as the research objects, and under the condition of disturbance ${\delta _2}$ and ${\delta _3}$ to varying degrees, the law of CMV infection was analyzed by numerical simulation. Selection of parameter $A = 2$, $\beta = 9.8 \times {10^{{\rm{ - }}3}}$, $\rho = 29$, $b = 0.6$ in Model (2), take $\alpha = 2.5$, ${b_1} = 0.1$ and $\alpha = 5, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {b_1} = 0.2$ separately. The corresponding basic regeneration numbers ${R_0} = 1.6917$ and $\hat R_0^s = 0.576$ can be calculated. At this time, the equilibrium point ${E^*}$and ${E_0}$ of the corresponding Model (1) are stable.

We can assume that the disturbance of the mortality rate of healthy cells remains unchanged, that is, ${\delta _1}$ is set as 0.02. And the disruption values of the mortality rates of infected cells and CMV are shown in Table 1. We can calculate that the threshold values $R_0^s$ is greater than 1 in all four cases. According to Theorem 2, CMV persists and vibrates around the equilibrium point. The amplitude of the vibration is still shown in table 1. Take initial value $T(0) = 1000$, $I(0) = 1000$, $V(0) = 1000$. The infected cells and CMV curves of Model (2) were simulated by MATLAB, as shown in Figures 1–4 below.

As shown in Figures 1–4 and Table 1, the solution of Model (2) vibrates around the solution of Model (1). Based on Figure 1, comparing Figures 2, 3 and 4, the amplitude of vibration increases with the increase of disturbance intensity.

Let $\alpha = 5$, ${b_1} = 0.2$, ${\delta _2} = 0.04$, ${\delta _3} = 0.04$, and the values of other parameters remain the same. We can calculate$\hat R_0^s < 1$. According to theorem 4, CMV will eventually become extinct. The initial value is selected as follows

The simulation of Model (2) is shown in Figure 5 below.

As can be seen from Figure 5, both infected cells and viruses are extinct, and normal cells fluctuate randomly.

It can be seen in the above simulation curves, the stability of the equilibrium point of Model (1) will not be changed by the perturbation of the parameters. The solution of Model (2) will vibrate nearby, and the magnitude of the vibration depends on the perturbation value of the parameters.

5.   Conclusion

In this paper, the sufficient conditions for the existence of asymptotic, ergodic and extinctive solutions were obtained by studying the dynamical properties of the cytomegalovirus model with random behavior. It was known that the ultimate trend of the solution is to oscillate around the deterministic equilibrium point, and the amplitude of the vibration is determined by the intensity of the disturbance. According to the clinical diagnosis standard of CMV infection, only by controlling the main factors, which are $\alpha$ and ${b_1}$ that cause the change of CMV quantity as much as possible, can CMV be exterminated or controlled in a controllable range, so as to ensure the health of children.

Acknowledgments

This work was supported by the project of Nature Scientific Foundation of Heilongjiang Province (A2016004), National Natural Science Foundation of China (11801122).

Conflict of interest

The authors declare there is no conflict of interest.