# Analysis of a Multiple Delays Model for Treatment of Cancer with Oncolytic Virotherapy

Authors: Adil El Alami laaroussi; Mohamed El Hia; Mostafa Rachik; Rachid Ghazzali
Journal: Computational and Mathematical Methods in Medicine
 (2019)

---

## Abstract

Despite advanced discoveries in cancerology, conventional treatments by surgery, chemotherapy, or radiotherapy remain ineffective in some situations. Oncolytic virotherapy, i.e., the involvement of replicative viruses targeting specific tumor cells, opens new perspectives for better management of this disease. Certain viruses naturally have a preferential tropism for the tumor cells; others are genetically modifiable to present such properties, as the lytic cycle virus, which is a process that represents a vital role in oncolytic virotherapy. In the present paper, we present a mathematical model for the dynamics of oncolytic virotherapy that incorporates multiple time delays representing the multiple time periods of a lytic cycle. We compute the basic reproductive ratioR0, and we show that there exist a disease-free equilibrium point (DFE) and an endemic equilibrium point (DEE). By formulating suitable Lyapunov function, we prove that the disease-free equilibrium (DFE) is globally asymptotically stable if R0<1 and unstable otherwise. We also demonstrate that under additional conditions, the endemic equilibrium is stable. Also, a Hopf bifurcation analysis of our dynamic system is used to understand how solutions and their stability change as system parameters change in the case of a positive delay. To illustrate the effectiveness of our theoretical results, we give numerical simulations for several scenarios.


---

## Body

## 1. Introduction

The continuous improvement of conventional treatments in cancerology (surgery, chemotherapy, and radiotherapy) allows for a major progress in the fight against cancer. Nevertheless, in some situations, these modes of therapy may be ineffective. The development of new therapeutic strategies, therefore, appears essential in order to improve the healing of this disease. Thus, in the last decades, virotherapy of cancers appears to be a credible alternative to some situations, due to advanced discoveries and accurate informations about viruses and also the production of recombinant viral vectors which can be used in cancer gene therapy. The use of replicative viruses as antitumor therapeutic agents (oncolytic viruses) is based on the idea that they reproduce preferentially within the tumor cells; however, normal cells remain immune to infection. These viruses (oncolytic viruses) are either virus with a natural ability to replicate preferentially within tumor cells or viruses genetically modified to hold this property. Genetic modifications are primarily used to improve the specificity of viruses against tumor cells, by targeting a particular surface molecule, by deleting specific viral genes required for replication in healthy cells, or by using activatable viral promoters only in tumor cells [1]. Using viruses to treat cancer is not a new concept. Viruses have attracted interest as anticancer therapeutics since the beginning of the 20th century. However, for several years, research in this field was limited due to technological limitations. In the last 30 years, by increasing understanding of the nature of viruses, their mechanisms of oncolytic activity and their ability to manipulate and exploit genetically has prompted a new wave of oncolytic virotherapy. Today, there is extensive literature describing progress in both theoretical and clinical trials of oncolytic viruses. For more details, we refer the interested reader to [2–7].Mathematical modeling of oncolytic virotherapy can illuminate the underlying dynamics of treatment systems and lead to optimal treatment strategies. Several studies have been the subject of the study of virotherapy. The first mathematical models of oncolytic viral therapy used ordinary differential equations to describe the fundamental interactions between two types of tumor cells (infected cells and uninfected cells) [8, 9]. Other works consider spatial representation of tumors [10, 11, 12], multiscale effects [13], and stochastic processes [14]. For a review of different mathematical modeling approaches ranging from ordinary differential equations to spatially explicit agent-based models, see [15–17].Biological experiments helped to understand and explain the lytic cycle, which takes place in six stages; the six stages are as follows: attachment, penetration, transcription, biosynthesis, maturation, and lysis. To infect a new cell, a virus must penetrate inside the cell through the plasma membrane; the virus attacks a receptor on the cell membrane and then releases its genetic material in the cell. In the third step, the host cell’s DNA is degraded and the cell’s metabolism is directed to initiate the fourth step; biosynthesis, here the virus uses cellular mechanisms to constitute a large amount of viral components and, in the meantime, destroys the DNA of the host cell. Then, it enters the last two stages, maturation and lysis. When many copies of viral are manufactured, they are assembled into complete formed viruses. About 25 minutes after initial infection, approximately 200 new bacteriophages (virions) are formed. Once enough virions have matured and accumulated, specialized viral proteins are used to dissolve the bacterial cell wall, where they can go on to infect other cells and another lytic cycle begins (for more details on the lytic cycle, see [5, 18]). In this work, the dynamics of oncolytic virotherapy are studied by incorporating the viral lytic cycle time. The duration of the intracellular viral life cycle is an essential factor in viral therapy. For example, some viruses require only 30 minutes, some viruses take several hours to complete this process, and some may take days [19]. Therefore, it is necessary and realistic to consider and taking into account the cycle time in modeling the oncolytic virotherapy which allows us to better predict its dynamics. We construct a mathematical model of virotherapy with multiple delays representing the six time periods of the lytic cycle; it is assumed that the time of each stage of the lytic cycle is constant.Several studies have studied and analyzed systems of delayed differential equations that model virotherapy. In the paper entitled “Hopf Bifurcation Analysis in a Delayed System for Cancer Virotherapy” [20], the authors consider a delayed differential equation system. In [19], Wang et al. propose a mathematical model for oncolytic virotherapy where they consider the time period of the viral lytic cycle as a delay parameter. The novelty of our work is modeling the variation of duration in the intracellular viral life cycle by adding multiple delays; each one represents the time period of each stage of the lytic cycle. We compute the basic reproductive ratio R0, and we show that there exist a disease-free equilibrium point (DFE) and an endemic equilibrium point (DEE). By formulating suitable Lyapunov function, we prove that the DFE is globally asymptotically stable if R0<1 and unstable otherwise. We also demonstrate that under additional conditions, the DEE is stable. Furthermore, a bifurcation analysis of our dynamical system is used to understand how solutions and their stability change as the parameters in the system vary. To illustrate our theoretical results, numerical simulations are also presented for several scenarios.This paper is organized as follows. In Section2, we present our mathematical model. In Section 3, we compute the equilibrium of our model and investigate its stability. Following that, in Section 4, a bifurcation analysis of the dynamical system is used to understand how the solutions and their stability change as the parameters change. Numerical studies are shown, in Section 5, to validate the analytical results. Finally, we conclude the paper in Section 6.


## 2. The Basic Mathematical Model

In a previous work [21], we analyzed the stability of a nonlinear system of differential equations based on the models proposed by [22, 23]. Our model contains three variables, which are, uninfected tumor cell population xt, infected tumor cell population yt, and free virus particles which are outside cell vt, and it has the following form:(1)dxtdt=rxt1−xt+ytK−βxtvt−ρxtyt,dytdt=βxtvt−δyt,dvtdt=bδyt−γvt−βxtvt,where x0=x0, y0=y0, and v0=v0 are given.The termrx1−x+y/K describes the logistic growth rate of an uninfected tumor cell population xt. The constant r>0 is the growth rate, with K being the carrying capacity or maximal tumor size so that x+y≤K. The term βxv represents the rate of infected cells by free virus vt, with β>0 being the corresponding constant rate. The term ρxy models infection from an encounter between an infected cell and an uninfected cell resulting in cell fusion that produces a syncytium, with ρ>0 being the constant rate describing cell to cell fusion with the formation of syncytia. Infected cells die at a rate of δy, and γv is the rate of elimination of free virus particles by various causes including nonspecific binding and generation of defective interfering particles. Its burst size models the virus replication ability, the burst size of a virus, which is an essential parameter of virus reproduction. So, our model includes also a parameter b that models the burst size.As we mentioned in Introduction, the model (1) did not take into account the time needed to complete the lytic cycle. As a reminder, the lytic cycle is the process where a virus overtakes a cell and uses the cellular machinery of its host to reproduce. Copies of the virus fill the cell to bursting, killing the cell and releasing viruses to infect more cells. The duration of this process varies from virus and more cells. Wang et al. [19] proposed a model of virotherapy with a single delay time; the originality of our work is to make a generalization by introducing 6 delays representing each period of stages of the lytic cycle in order to describe a more realistic situation because the virus goes through 6 stages of life and each one of them may have a different delay. We denote τii=1,2,3,4,5,n=6, the different times period of the lytic cycle. The rate of change of infected tumor cells at time t will be determined by the tumor cell population and free virus at time t−τi, namely, xt−τivt−τi; for more details about the lytic cycle, see Figure 1. Therefore, the model we propose is given as follows:(2)dxtdt=rxt1−xt+ytK−βxtvt−ρxtyt,dytdt=∑i=1nβixt–τivt–τi−δyt,dvtdt=bδyt−γvt−βxtvt,where τ0=0, τ1<τ2<⋯<τn, and β=∑i=1nβi.Figure 1
Schematic diagram of the model for virotherapy.Using the Van Den Driesseche and Watmough next-generation approach, we calculate the basic reproductive ratio of system (2), which leads to(3)R0=βKbβK+γ.This parameter plays a major role in our analysis. It represents the number of new virus particles generated by a single virus particle that is inserted into a tumor consisting entirely of uninfected tumor cells [12].


## 3. Model Analysis

In this section, we show the existence of the equilibrium points and we study their stabilities. System (2) has three equilibrium, E0=0,0,0, E1=K,0,0, and the positive equilibrium E∗x∗,y∗,v∗, where(4)x∗=γβb−1,y∗=γrKbβ−Kβ−γβb−1rγ+δKβb−1+ρKγ,v∗=δrKbβ−Kβ−γβrγ+δKβb−1+ρKγ.The Jacobian matrix of system (2) at an arbitrary point is given by(5)J=−rx∗K−ρ+rKx∗βx∗∑i=1nβiv∗e−λτi−δ∑i=1nβix∗e−λτi−βv∗bδ−γ−βx∗.The first equilibrium+ represents the total success of therapy. It is easy to prove that E0 is always unstable. Biologically, the instability of this equilibrium is because, in the absence of the virus, the number of infected cells y will remain at 0, and tumor cell population increases. The second equilibrium E1 represents the failure of virotherapy, as the tumor achieved its maximal size K. The partial success of virotherapy is represented by the third equilibrium E∗. The approach that we use to prove the stability of the steady states is divided into two parts: the first one concerns the necessary condition of stability when there is no delay τi=0 for i=1,…,n. In the second step, we prove that the matrix (5) does not have any imaginary eigenvalue. However, in our case, it was not easy to apply the classical theorems of stability because we deal with a system with multiple discrete delays as a summation ∑i=1nβixt−τivt−τi. To solve this problem, we brought the lemma below which allowed us to write the characteristic equation of (5) in a suitable form that allows the application of classical stability results.Lemma 1.
Forai∈R, we have(6)∑i=1nai2=∑i=1nai2+2∑i=0n−1∑j=i+1naiaj.Proof.
This result can be proved easily by induction.Remark 1.
We note that our approach has the limit to be specific to the model (2); it may not be appropriate for other models. The extension of our method to other models could be considered as one of the perspectives of this work.


### 3.1. Free Equilibrium

System (2) always has a disease-free equilibrium in the form E1=K,0,0.Proposition 1.
IfR0<1, then E1 is locally asymptotically stable.Proof.
The Jacobian matrix evaluated atE1 is(7)JE1=−r−ρK−r−βK0−δ∑i=1nβiKe−λτi0bδ−γ−βK.
The characteristic polynomial ofJE1 is(8)PE1λ=−r+λλ2+λγ+βK+δ+δβK+γ−∑i=1nβiKbδe−λτi.
Ifτ1=τ2=⋯=τn=0, then(9)PE1λ=−r+λλ2+λγ+βK+δ+δβK+γ−βKbδ.
In this case, the eigenvalues of the matrixJE1 are(10)λ1=−r,λ2=−γ+βK+δ−Δ2,λ3=−γ+βK+δ+Δ2,where Δ=γ+βK−δ2+4βKbδ.
The eigenvaluesλ1 and λ2 are both negatives for all non-negative parameter values, while the eigenvalue λ3 can be negative, positive, and zero. For R0<1, we have(11)R0<1⟺βKbβK+γ<1⟺4βKbδ<4δβK+γ⟺4βKbδ+γ+βK–δ2<4δβK+γ+γ+βK–δ2⟺Δ<γ+βK+δ2⟺λ3=–γ+βK+δ+Δ2<0.
Hence, all three eigenvalues are negatives. SoE1 is locally asymptotically stable when τi=0 for i=1,…,n.
Now, ifτii=1,…,n are arbitrary and as λ=−r is a root of equation (8), we only need to consider(12)λ2+λγ+βK+δ+δβK+γ−∑i=1nβiKbδe−λτi=0,which is equivalent to(13)λ2+λγ+βK+δ+δβK+γ=∑i=1nβiKbδe−λτi.
Ifλ=ωi is a root of equation (12), after substituting and separating real and imaginary parts, we have(14)−ω2+δβK+γ=Kbδ∑i=1nβicosωτi,−ωγ+βK+δ=Kbδ∑i=1nβisinωτi.
Adding the squares of both equations from (14), one has(15)ω4+δ+βK+γ2−2δβK+γω2+δ2γ+βK2=Kbδ2∑i=1nβicosωτi2+∑i=1nβisinωτi2.
Using Lemma 1 and after algebraic manipulations, equation (15) can also be written in the following form:(16)ω4+a1ω2+a2+a3=0,where(17)a1=γ+βK2+δ2,a2=2bKδ2∑i=0n−1∑j=i+1nβiβj1−cosωτi−τj,a3=δ2γ+βK2−βbKδ2=δ2γ+βK21−R02.
We havea1>0 and a2>0, and when R0<1, a3>0. Therefore, there is no root λ=ωi, with ω≥0 or equation (12), implying that the roots of equation (12) cannot cross the purely imaginary axis. Thus, all roots of equation (12) have a negative part. Then, the equilibrium point E1 is locally asymptotically stable.
By using a Lyapunov function, we will prove that the equilibrium pointE1 is globally asymptotically stable when R0<1. To study the dynamics of system (2) when τi≥0i=1,…,n, we need to consider a suitable phase space. For τn>0, we denote by C=C−τn;0;R3 the Banach space of continuous functions mapping the interval −τn;0 into R3 with the norm φθ=sup−rn≤θ≤0φθ for φ∈C. The non-negative cone of C is denoted by C+=C−τn,0,R+3.Theorem 1.
IfR0<1, then E1 is globally asymptotically stable.Proof.
Letφ=φ1,φ2,φ3=x,y,v with xθ=φ1θ,yθ=φ2θ,vθ=φ3θ for θ∈−τn,0, consider a Lyapunov function given by(18)Vφ=bφ20+φ30+b∑i=1nβi∫−τi0φ1sφ3sds.
The derivative along a solution is given by(19)V˙φ=b∑i=1nβiφ1−τiφ3−τi−δφ20+bδφ20−βφ10φ30−γφ30+b∑i=1nβiφ10φ30−φ1−τiφ3−τi=bβφ10φ30−βφ10φ30−γφ30≤bβ–βK–γφ30≤bβK1–1R0φ30,when R0<1, we have V˙φ≤0. If V˙φ=0, then φ∈R3/V˙φ=0=E1. The classical LaSalle’s invariance Principle implies that E1 is globally attractive. This confirms the globally asymptotical stability of E1.


### 3.2. Endemic Equilibrium

Here, we study the stability of the endemic equilibrium pointE∗.Theorem 2.
Equilibrium pointE∗ is locally asymptotically stable for τi≥0i=1,…,n if the following assumptions are satisfied:(20)A1R0>1,A2fb<rγ+Kργ–βKδ,A3Ki>0,fori=1,…,3,where(21)fb=rγ+δKβb−β+ρKγγKβb−1Kβδb−1+Kβγ+rγ−δb–1+γbδb–1Kβb–δ–βK–γ,K1=A2−2B˜+β0D,K2=B˜+β0D2−2AC+β0E−∑i=1nβiD2+2∑i=1n−1∑j=i+1nβiβjD21−cosωτi−τj,K3=C+β0E2−∑i=1nβiE2+2∑i=1n−1∑j=i+1nβiβjE21–cosωτi−τj,and(22)A=βx∗+γ+δ+rKx∗,B˜=rKx∗βx∗+γ+δ+δβx∗+γ−δβ2x∗v∗,C=rKx∗δβx∗+γ−δβ2x∗v∗,D=rK+ρx∗v∗−bδx∗,E=rK+ργ+βδbx∗v∗−rKbδx∗2.Proof.
The Jacobian matrix atE∗ is given by(23)JE∗=−rKx∗−ρ+rKx∗−βx∗∑i=1nβiv∗e−λτi−δ∑i=1nβiv∗e−λτi−βv∗bδ−βx∗−γ.
The characteristic equation associated withJE∗ is given by(24)λ3+Aλ2+B˜λ+C+∑i=1nβie−λτiDλ+E=0,where A,B˜,C,D, and E are defined as in Theorem 2.
Consideringτ1=τ2=⋯=τn=0, equation (23) becomes(25)JE1=−rKγβbγ−1−ρ+rKγβbγ−1−γb−1βv∗−δγb−1−βv∗bδ−γbb−1,and equation (24) becomes(26)λ3+b1λ2+b2λ+b3=0,with(27)b1=δ+γβKb−1+βγK+rγKβb–1,b2=rγδb−δ+γb2+γbβKb–12+δrKβb–βK–γrγβ+βδKβb–β+pKγ×rγ+Kργ−βγKKb−1,b3=δrγKβb–βK–γKβb–1.
By the Routh–Hurwitz Criterion, all roots of the polynomial (26) have negative real parts if and only if H1=b1>0, H2=b1b2−b3>0, and H3=b2H2>0. When R0>1, we have H1=b1>0 and b2>0. Since H3=b3H2, we only need to consider H2=b1b2−b3.
After some algebraic manipulations ([21]) we can prove that(28)H2>0⟺fb<rγ+Kργ−βKδ.
So we conclude that whenR0>1 and fb<rγ+Kργ−βKδ, the endemic equilibrium is locally asymptotically stable for τ1=⋯=τn=0.
Consider now the case whenτ1,…,τn are arbitrary. Finding roots of the equation (24) is impossible explicitly. Instead, we look for the condition under which it has no purely imaginary roots.
Letλ=ωiω>0 be a purely imaginary roots of (24), then(29)−ω3i−Aω2+B˜ωi+C+∑i=1nβiDωi+Ee−iωτi=0,which is equivalent to(30)−ω3i−Aω2+B˜ωi+C+∑i=1nβiDωi+Ecosωτi−isinωτi=0.
Separating real and imaginary parts leads to(31)Aω2−C−β0E=∑i=1nβiEcosωτi+Dωsinωτi,ω3−B˜ω−β0Dω=∑i=1nβiDωcosωτi−Esinωτi.
Adding the squares of both equations together gives(32)ω6+A2–2B˜+β0Dω4+B˜+β0D2–2AC+β0Eω2+C+β0E2=–∑i=1nβiEcosωτi+Dωsinωτi2=–∑i=1nβiEcosωτi+Dωsinωτi2+–∑i=1nβiDωcosωτi−Esinωτi2.
Using Lemma 1 and after some algebraic manipulations, equation (32) can also be written in the following form:(33)ω6+K1ω4+K2ω2+K3=0,where Kii=1,2,3 are as in Theorem 2. If Ki>0, then all roots of (24) have negative real parts. Hence, the proof is complete.


## 3.1. Free Equilibrium

System (2) always has a disease-free equilibrium in the form E1=K,0,0.Proposition 1.
IfR0<1, then E1 is locally asymptotically stable.Proof.
The Jacobian matrix evaluated atE1 is(7)JE1=−r−ρK−r−βK0−δ∑i=1nβiKe−λτi0bδ−γ−βK.
The characteristic polynomial ofJE1 is(8)PE1λ=−r+λλ2+λγ+βK+δ+δβK+γ−∑i=1nβiKbδe−λτi.
Ifτ1=τ2=⋯=τn=0, then(9)PE1λ=−r+λλ2+λγ+βK+δ+δβK+γ−βKbδ.
In this case, the eigenvalues of the matrixJE1 are(10)λ1=−r,λ2=−γ+βK+δ−Δ2,λ3=−γ+βK+δ+Δ2,where Δ=γ+βK−δ2+4βKbδ.
The eigenvaluesλ1 and λ2 are both negatives for all non-negative parameter values, while the eigenvalue λ3 can be negative, positive, and zero. For R0<1, we have(11)R0<1⟺βKbβK+γ<1⟺4βKbδ<4δβK+γ⟺4βKbδ+γ+βK–δ2<4δβK+γ+γ+βK–δ2⟺Δ<γ+βK+δ2⟺λ3=–γ+βK+δ+Δ2<0.
Hence, all three eigenvalues are negatives. SoE1 is locally asymptotically stable when τi=0 for i=1,…,n.
Now, ifτii=1,…,n are arbitrary and as λ=−r is a root of equation (8), we only need to consider(12)λ2+λγ+βK+δ+δβK+γ−∑i=1nβiKbδe−λτi=0,which is equivalent to(13)λ2+λγ+βK+δ+δβK+γ=∑i=1nβiKbδe−λτi.
Ifλ=ωi is a root of equation (12), after substituting and separating real and imaginary parts, we have(14)−ω2+δβK+γ=Kbδ∑i=1nβicosωτi,−ωγ+βK+δ=Kbδ∑i=1nβisinωτi.
Adding the squares of both equations from (14), one has(15)ω4+δ+βK+γ2−2δβK+γω2+δ2γ+βK2=Kbδ2∑i=1nβicosωτi2+∑i=1nβisinωτi2.
Using Lemma 1 and after algebraic manipulations, equation (15) can also be written in the following form:(16)ω4+a1ω2+a2+a3=0,where(17)a1=γ+βK2+δ2,a2=2bKδ2∑i=0n−1∑j=i+1nβiβj1−cosωτi−τj,a3=δ2γ+βK2−βbKδ2=δ2γ+βK21−R02.
We havea1>0 and a2>0, and when R0<1, a3>0. Therefore, there is no root λ=ωi, with ω≥0 or equation (12), implying that the roots of equation (12) cannot cross the purely imaginary axis. Thus, all roots of equation (12) have a negative part. Then, the equilibrium point E1 is locally asymptotically stable.
By using a Lyapunov function, we will prove that the equilibrium pointE1 is globally asymptotically stable when R0<1. To study the dynamics of system (2) when τi≥0i=1,…,n, we need to consider a suitable phase space. For τn>0, we denote by C=C−τn;0;R3 the Banach space of continuous functions mapping the interval −τn;0 into R3 with the norm φθ=sup−rn≤θ≤0φθ for φ∈C. The non-negative cone of C is denoted by C+=C−τn,0,R+3.Theorem 1.
IfR0<1, then E1 is globally asymptotically stable.Proof.
Letφ=φ1,φ2,φ3=x,y,v with xθ=φ1θ,yθ=φ2θ,vθ=φ3θ for θ∈−τn,0, consider a Lyapunov function given by(18)Vφ=bφ20+φ30+b∑i=1nβi∫−τi0φ1sφ3sds.
The derivative along a solution is given by(19)V˙φ=b∑i=1nβiφ1−τiφ3−τi−δφ20+bδφ20−βφ10φ30−γφ30+b∑i=1nβiφ10φ30−φ1−τiφ3−τi=bβφ10φ30−βφ10φ30−γφ30≤bβ–βK–γφ30≤bβK1–1R0φ30,when R0<1, we have V˙φ≤0. If V˙φ=0, then φ∈R3/V˙φ=0=E1. The classical LaSalle’s invariance Principle implies that E1 is globally attractive. This confirms the globally asymptotical stability of E1.


## 3.2. Endemic Equilibrium

Here, we study the stability of the endemic equilibrium pointE∗.Theorem 2.
Equilibrium pointE∗ is locally asymptotically stable for τi≥0i=1,…,n if the following assumptions are satisfied:(20)A1R0>1,A2fb<rγ+Kργ–βKδ,A3Ki>0,fori=1,…,3,where(21)fb=rγ+δKβb−β+ρKγγKβb−1Kβδb−1+Kβγ+rγ−δb–1+γbδb–1Kβb–δ–βK–γ,K1=A2−2B˜+β0D,K2=B˜+β0D2−2AC+β0E−∑i=1nβiD2+2∑i=1n−1∑j=i+1nβiβjD21−cosωτi−τj,K3=C+β0E2−∑i=1nβiE2+2∑i=1n−1∑j=i+1nβiβjE21–cosωτi−τj,and(22)A=βx∗+γ+δ+rKx∗,B˜=rKx∗βx∗+γ+δ+δβx∗+γ−δβ2x∗v∗,C=rKx∗δβx∗+γ−δβ2x∗v∗,D=rK+ρx∗v∗−bδx∗,E=rK+ργ+βδbx∗v∗−rKbδx∗2.Proof.
The Jacobian matrix atE∗ is given by(23)JE∗=−rKx∗−ρ+rKx∗−βx∗∑i=1nβiv∗e−λτi−δ∑i=1nβiv∗e−λτi−βv∗bδ−βx∗−γ.
The characteristic equation associated withJE∗ is given by(24)λ3+Aλ2+B˜λ+C+∑i=1nβie−λτiDλ+E=0,where A,B˜,C,D, and E are defined as in Theorem 2.
Consideringτ1=τ2=⋯=τn=0, equation (23) becomes(25)JE1=−rKγβbγ−1−ρ+rKγβbγ−1−γb−1βv∗−δγb−1−βv∗bδ−γbb−1,and equation (24) becomes(26)λ3+b1λ2+b2λ+b3=0,with(27)b1=δ+γβKb−1+βγK+rγKβb–1,b2=rγδb−δ+γb2+γbβKb–12+δrKβb–βK–γrγβ+βδKβb–β+pKγ×rγ+Kργ−βγKKb−1,b3=δrγKβb–βK–γKβb–1.
By the Routh–Hurwitz Criterion, all roots of the polynomial (26) have negative real parts if and only if H1=b1>0, H2=b1b2−b3>0, and H3=b2H2>0. When R0>1, we have H1=b1>0 and b2>0. Since H3=b3H2, we only need to consider H2=b1b2−b3.
After some algebraic manipulations ([21]) we can prove that(28)H2>0⟺fb<rγ+Kργ−βKδ.
So we conclude that whenR0>1 and fb<rγ+Kργ−βKδ, the endemic equilibrium is locally asymptotically stable for τ1=⋯=τn=0.
Consider now the case whenτ1,…,τn are arbitrary. Finding roots of the equation (24) is impossible explicitly. Instead, we look for the condition under which it has no purely imaginary roots.
Letλ=ωiω>0 be a purely imaginary roots of (24), then(29)−ω3i−Aω2+B˜ωi+C+∑i=1nβiDωi+Ee−iωτi=0,which is equivalent to(30)−ω3i−Aω2+B˜ωi+C+∑i=1nβiDωi+Ecosωτi−isinωτi=0.
Separating real and imaginary parts leads to(31)Aω2−C−β0E=∑i=1nβiEcosωτi+Dωsinωτi,ω3−B˜ω−β0Dω=∑i=1nβiDωcosωτi−Esinωτi.
Adding the squares of both equations together gives(32)ω6+A2–2B˜+β0Dω4+B˜+β0D2–2AC+β0Eω2+C+β0E2=–∑i=1nβiEcosωτi+Dωsinωτi2=–∑i=1nβiEcosωτi+Dωsinωτi2+–∑i=1nβiDωcosωτi−Esinωτi2.
Using Lemma 1 and after some algebraic manipulations, equation (32) can also be written in the following form:(33)ω6+K1ω4+K2ω2+K3=0,where Kii=1,2,3 are as in Theorem 2. If Ki>0, then all roots of (24) have negative real parts. Hence, the proof is complete.


## 4. Hopf Bifurcation

In this section, we will study the Hopf bifurcation of system (2) but only in the case of one positive term of delay. In fact, it is too difficult to study the general case with n>1, which can be considered as a perspective of this work. Consider n=1, then (24) becomes(34)λ3+Aλ2+B˜λ+C+β0Dλ+E=−β1Dλ+Ee–λτ1,if λ=ωi is a root of (34). After substituting and separating real and imaginary parts, we have(35)−Aω2+C+β0E=−β1Ecosωτ1+Dωsinωτ1,−ω3+B˜ω+β0Dω=−β1Dωcosωτ1−Esinωτ1.Adding the squares of both equations together gives(36)ω6+K1ω4+K2ω2+K3=0,where K1, K2, and K3 are as follows:(37)K1=A2−2B˜+β0D,K2=B˜+β0D2−2AC+β0E−β1D2,K3=C+β0E2−β1E2.Denoteω0 the biggest positive root of (37); then from (35), we have(38)cosω0τ1=Dω04+AE–DB˜+β0Dω02–C+β0EEβ1Dω02+E2.Let(39)τ1,0j=1ω0arcosDω04+AE–DB˜+β0Dω02–C+β0EEβ1Dω02+E2+2jπ,j=1,2,….Then, we can defineτ∗=minj≥1τ1,0j as the first value of τ1 when characteristic roots cross the imaginary axis.Further, differentiating equation (34) with respect to τ1, we get(40)3λ2+2Aλ+B˜+β0D+β1De−λτ1−β1Dλ+Ee−λτ1dλdτ1=β1Dλ+Ee−λτ1.This gives(41)dλdτ1−1=3λ2+2Aλ+B˜+β0Dβ1Dλ+Ee−λτ1+DλDλ+E−τ1λ,and after some algebraic manipulations, we get(42)dλdτ1−1=2λ3+Aλ2−C−β0E−λ2λ3+Aλ2+B˜λ+C+β0Dλ+E+Eλ2Dλ+E−τ1λ.Thus,(43)dλdτ1−1λ=ω0i=−i2ω03−Aω02−C−β0Eω02−iω03−Aω02+B˜ω0i+C+β0Dω0i+E+−Eω02Dω0i+E−τ1ω0i=−Aω02+C+β0E−i2ω03ω02−Aω02+C+β0E+i−ω03+B˜+βDω0+−Eω02E+Dω0i−τ1ω0i,and(44)Redλdτ1−1λ=ω0i=2ω06−A2−2B˜+β0Dω04−β12E2−C+β0E2ω02−Aω02+C+β0E2+−ω03+B˜+β0Dω02.So, ifn=1 and K3<0, there exists a Hopf bifurcation as Redλ/dτ1−1λ=ω0i>0.In conclusion, we have the following Hopf bifurcation result.Theorem 3.
In the case where system (2) has only one no zero delays and if K3<0, then system (2) undergoes a Hopf bifurcation at the endemic equilibrium.


## 5. Numerical Results

In this section, we present numerical simulations to illustrate the various theoretical results previously obtained. Thus, we draw first the curves of system (2) for parameters verifying R0 less than 1, and we shall do the same for parameters verifying R0 upper to 1. All simulations are performed using the parameter values in Table 1 are taken from [22].Table 1
Model parameters.
Height parameters
Descriptions
Values
r
Growth rate constant
0.206
K
Maximal tumor size
2139
β
i
Infection rate
Variables
ρ
Cell to cell fusion rate constant
0.2145
δ
Infected cells death rate
0.5115
b
Burst size of a virus
Variables
γ
Elimination rate of free virus particles
0.001
τ
i
Time to complete lytic cycle variables
VariablesSince our model considers population of cells, we convert tumor volume to cell population by assuming1mm3 corresponds to 106 cells [22]. For our numerical simulation, we consider cell populations x and y, virus population v is expressed in units of 106, and using the same manner as in [22], we assume that the tumor is completely eliminated, which indicates the total success of the virotherapy, when the total population of tumor cells is reduced to one cell, which means that in the adopted units ut=xt+yt=106.Figure2 presents the curves of system (2), using various initial conditions when n=4 and R0=0.7. In Section 3, using a suitable Lyapunov function, we have proved that, in this case, R0<1, the disease-free equilibrium E1 is globally asymptotically stable. From this figure, we see that the curves converge to the free equilibrium E1, that is the virotherapy fails as the population of tumor cells increase and the population of infected tumor decrease.Figure 2
The lytic cycle of oncolytic viruses.Figure3 provides the curves of system (2) using various initial conditions when n=4, R0=1.5, and the other conditions of Theorem 2 are satisfied. We have theoretically proved in Section 3 by using the technique of stability in a delayed system that the endemic equilibrium E∗ is locally asymptotically stable. From this figure, we see that the curves converge to positive and finite limit, which is the endemic equilibrium. The stability of the equilibrium E∗ implies that a permanent reduction of the tumor load can be reached, even if the virotherapy does not succeed completely.Figure 3
Dynamics of virotherapy whenR0=0.7, β0=10−3, β1=10−5, β2=2×10−5, β3=3×10−5, β4=4×10−5, τ1=0.2, τ2=1, τ3=2, and τ4=3. The initial conditions are x0=127, y0=0, and v0=30.Figures4–6 show that for τ<τ∗, the equilibrium point E∗ is asymptotically stable and, for τ>τ∗, the equilibrium point E∗ is unstable, and when τ=τ∗, a Hopf bifurcation of periodic solutions of system (2) occurs at E∗ (Figure 7).Figure 4
Dynamics of virotherapy whenR0=1,5, β0=10−3, β1=10−5, β2=2×10−5, β3=3×10−5, β4=4×10−5, τ1=0.2, τ2=1, τ3=2, and τ4=3. The initial conditions are x0=127, y0=0, and v0=30.Figure 5
Dynamics of virotherapy whenR0=1,5, β0=10−5, β1=10−4, and τ=7<τ∗=8.8368. The initial conditions are x0=127, y0=0, and v0=30.Figure 6
Dynamics of virotherapy whenR0=1,5, β0=10−5, β1=10−4, and τ=τ∗=8.8368. The initial conditions are x0=127, y0=0, and v0=30.Figure 7
Dynamics of virotherapy whenR0=1,5, β0=10−5, β1=10−4, and τ=10>τ∗=8.8368. The initial conditions are x0=127, y0=0, and v0=30.


## 6. Conclusion

The work in this paper contributes to a growing literature on modeling oncolytic virotherapy; we present a mathematical model for the dynamic of oncolytic virotherapy that incorporates multiple time delays representing the multiple time periods to complete the lytic cycle. We give the basic reproductive ratioR0, and we use it to investigate the stability of the equilibrium states. We prove by formulating suitable Lyapunov function that the disease-free equilibrium is globally asymptotically stable if the basic infection reproduction number R0<1, and when R0>1, the local stability of the endemic equilibrium point depends on function fb, representing the replication of the virus in virotherapy and other conditions. Furthermore, we show that there exists a bifurcation value for the lytic cycle period τ∗. For this, if τ<τ∗, the positive equilibrium endemic is locally asymptotically stable. The system undergoes a Hopf bifurcation around τ=τ∗ and when τ>τ∗, the system is unstable. The numerical simulation provides that if R0<1, the virotherapy fails as the population of tumors cells increases and the population of infected tumor decreases, and if R0>1, the virotherapy success and treatment will reach the equilibrium point endemic. The approach that we have introduced with multiple delays is specific to our model or to similar models in other fields. The incorporation of delay from a system that describes virotherapy is an interesting and realistic strategy, and several studies have adopted this method, for example, Wang’s work [24].



---

*Source: 1732815-2019-09-30.xml*