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

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 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 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.


Introduction
e 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. e development of new therapeutic strategies, therefore, appears essential in order to improve the healing of this disease. us, 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. e 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. ese 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][3][4][5][6][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. e 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][16][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. en, 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. e 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]. erefore, 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. e 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 R 0 , 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 R 0 < 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.
is paper is organized as follows. In Section 2, 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.

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 x(t), infected tumor cell population y(t), and free virus particles which are outside cell v(t), and it has the following form: where x(0) � x 0 , y(0) � y 0 , and v(0) � v 0 are given. e term rx(1 − ((x + y)/K)) describes the logistic growth rate of an uninfected tumor cell population x(t). e constant r > 0 is the growth rate, with K being the carrying capacity or maximal tumor size so that x + y ≤ K. e term βxv represents the rate of infected cells by free virus v(t), with β > 0 being the corresponding constant rate. e 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 cv 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. 2 Computational and Mathematical Methods in Medicine 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. e 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 τ i (i � 1, 2, 3, 4, 5, n � 6), the different times period of the lytic cycle. e 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, for more details about the lytic cycle, see Figure 1. erefore, the model we propose is given as follows: where τ 0 � 0, τ 1 < τ 2 < · · · < τ n , and β � n i�1 β i . Using the Van Den Driesseche and Watmough nextgeneration approach, we calculate the basic reproductive ratio of system (2), which leads to is 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].

Model Analysis
In this section, we show the existence of the equilibrium points and we study their stabilities. System (2) has three equilibrium, E 0 � (0, 0, 0), E 1 � (K, 0, 0), and the positive equilibrium , e Jacobian matrix of system (2) at an arbitrary point is given by e first equilibrium + represents the total success of therapy. It is easy to prove that E 0 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. e second equilibrium E 1 represents the failure of virotherapy, as the tumor achieved its maximal size K. e partial success of virotherapy is represented by the third equilibrium E * . e 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 n 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.
Proof. is 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. e extension of our method to other models could be considered as one of the perspectives of this work.
e Jacobian matrix evaluated at E 1 is Computational and Mathematical Methods in Medicine 3 e characteristic polynomial of J E 1 is If τ 1 � τ 2 � · · · � τ n � 0, then In this case, the eigenvalues of the matrix J E 1 are where Δ � (c + βK − δ) 2 + 4βKbδ. e eigenvalues λ 1 and λ 2 are both negatives for all nonnegative parameter values, while the eigenvalue λ 3 can be negative, positive, and zero. For R 0 < 1, we have Hence, all three eigenvalues are negatives. So E 1 is locally asymptotically stable when τ i � 0 for (i � 1, . . . , n). Now, if τ i (i � 1, . . . , n) are arbitrary and as λ � − r is a root of equation (8), we only need to consider which is equivalent to If λ � ωi is a root of equation (12), after substituting and separating real and imaginary parts, we have x

Computational and Mathematical Methods in Medicine
Adding the squares of both equations from (14), one has Using Lemma 1 and after algebraic manipulations, equation (15) can also be written in the following form: where We have a 1 > 0 and a 2 > 0, and when R 0 < 1, a 3 > 0. erefore, there is no root λ � ωi, with ω ≥ 0 or equation (12), implying that the roots of equation (12) cannot cross the purely imaginary axis. us, all roots of equation (12) have a negative part. en, the equilibrium point E 1 is locally asymptotically stable.
By using a Lyapunov function, we will prove that the equilibrium point E 1 is globally asymptotically stable when R 0 < 1. To study the dynamics of system (2) when τ i ≥ 0(i � 1, . . . , n), we need to consider a suitable phase space. For τ n > 0, we denote by C � C([− τ n ; 0]; R 3 ) the Banach space of continuous functions mapping the interval [− τ n ; 0] into R 3 with the norm ||φ(θ)|| � sup − r n ≤ θ ≤ 0 | φ(θ)| for φ ∈ C. e non-negative cone of C is denoted by , consider a Lyapunov function given by e derivative along a solution is given by e classical LaSalle's invariance Principle implies that E 1 is globally attractive. is confirms the globally asymptotical stability of E 1 .

Endemic Equilibrium.
Here, we study the stability of the endemic equilibrium point E * .

Theorem 2.
Equilibrium point E * is locally asymptotically stable for τ i ≥ 0(i � 1, . . . , n) if the following assumptions are satisfied: where

Computational and Mathematical Methods in Medicine
and Proof. e Jacobian matrix at E * is given by e characteristic equation associated with J E * is given by where A, B, C, D, and E are defined as in eorem 2.
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 which is equivalent to Separating real and imaginary parts leads to  Adding the squares of both equations together gives Using Lemma 1 and after some algebraic manipulations, equation (32) can also be written in the following form: where K i (i � 1, 2, 3) are as in eorem 2. If K i > 0, then all roots of (24) have negative real parts. Hence, the proof is complete.

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 if λ � ωi is a root of (34). After substituting and separating real and imaginary parts, we have Adding the squares of both equations together gives where K 1 , K 2 , and K 3 are as follows: Denote ω 0 the biggest positive root of (37); then from (35), we have en, we can define τ * � min j≥1 τ j 1,0 as the first value of τ 1 when characteristic roots cross the imaginary axis.
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 K 3 < 0, then system (2) undergoes a Hopf bifurcation at the endemic equilibrium.

Numerical Results
In this section, we present numerical simulations to illustrate the various theoretical results previously obtained.
us, we draw first the curves of system (2) for parameters verifying R 0 less than 1, and we shall do the same for parameters verifying R 0 upper to 1. All simulations are performed using the parameter values in Table 1 are taken from [22].
Since our model considers population of cells, we convert tumor volume to cell population by assuming 1 mm 3 corresponds to 10 6 cells [22]. For our numerical simulation, we consider cell populations x and y, virus population v is expressed in units of 10 6 , 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 u(t) � x(t) + y(t) � 10 6 . Figure 2 presents the curves of system (2), using various initial conditions when n � 4 and R 0 � 0.7. In Section 3, using a suitable Lyapunov function, we have proved that, in this case, R 0 < 1, the disease-free equilibrium E 1 is globally asymptotically stable. From this figure, we see that the curves converge to the free equilibrium E 1 , that is the virotherapy fails as the population of tumor cells increase and the population of infected tumor decrease. Figure 3 provides the curves of system (2) using various initial conditions when n � 4, R 0 � 1.5, and the other conditions of eorem 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. e 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. Figures 4-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).

Conclusion
e 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 ratio R 0 , 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 R 0 < 1, and when R 0 > 1, the local stability of the endemic equilibrium point depends on function f(b), 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. e system undergoes a Hopf bifurcation around τ � τ * and when τ > τ * , the system is unstable. e numerical simulation provides that if R 0 < 1, the virotherapy fails as the population of tumors cells increases and the population of infected tumor decreases, and if R 0 > 1, the virotherapy success and treatment will reach the equilibrium point endemic. e approach that we have introduced with multiple delays is specific to our model or to similar models in other fields. e 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].

Data Availability
Te disciplinary data used to support the findings of this study have been deposited in the Network Repository (http://www. networkrepository.com). Computational and Mathematical Methods in Medicine 11