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 ratio , 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 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.
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 , 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 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 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.
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 , infected tumor cell population , and free virus particles which are outside cell , and it has the following form:where , , and are given.
The term describes the logistic growth rate of an uninfected tumor cell population . The constant is the growth rate, with K being the carrying capacity or maximal tumor size so that . The term represents the rate of infected cells by free virus , with being the corresponding constant rate. The term models infection from an encounter between an infected cell and an uninfected cell resulting in cell fusion that produces a syncytium, with being the constant rate describing cell to cell fusion with the formation of syncytia. Infected cells die at a rate of , and 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 , 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 , namely, ; for more details about the lytic cycle, see Figure 1. Therefore, the model we propose is given as follows:where , , and .
Using the Van Den Driesseche and Watmough next-generation approach, we calculate the basic reproductive ratio of system (2), which leads to
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, , , and the positive equilibrium , where
The Jacobian matrix of system (2) at an arbitrary point is given by
The first equilibrium represents the total success of therapy. It is easy to prove that 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 represents the failure of virotherapy, as the tumor achieved its maximal size K. The partial success of virotherapy is represented by the third equilibrium . 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 for . 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 . 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. For , we have
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 .
Proposition 1. If , then is locally asymptotically stable.
Proof. The Jacobian matrix evaluated at isThe characteristic polynomial of isIf , thenIn this case, the eigenvalues of the matrix arewhere .
The eigenvalues and are both negatives for all non-negative parameter values, while the eigenvalue can be negative, positive, and zero. For , we haveHence, all three eigenvalues are negatives. So is locally asymptotically stable when for .
Now, if are arbitrary and as is a root of equation (8), we only need to considerwhich is equivalent toIf is a root of equation (12), after substituting and separating real and imaginary parts, we haveAdding the squares of both equations from (14), one hasUsing Lemma 1 and after algebraic manipulations, equation (15) can also be written in the following form:whereWe have and , and when , . Therefore, there is no root , with 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 is locally asymptotically stable.
By using a Lyapunov function, we will prove that the equilibrium point is globally asymptotically stable when . To study the dynamics of system (2) when , we need to consider a suitable phase space. For , we denote by the Banach space of continuous functions mapping the interval into with the norm for . The non-negative cone of C is denoted by .
Theorem 1. If , then is globally asymptotically stable.
Proof. Let with for , consider a Lyapunov function given byThe derivative along a solution is given bywhen , we have . If , then . The classical LaSalle’s invariance Principle implies that is globally attractive. This confirms the globally asymptotical stability of .
3.2. Endemic Equilibrium
Here, we study the stability of the endemic equilibrium point .
Theorem 2. Equilibrium point is locally asymptotically stable for if the following assumptions are satisfied:whereand
Proof. The Jacobian matrix at is given byThe characteristic equation associated with is given bywhere , and E are defined as in Theorem 2.
Considering , equation (23) becomesand equation (24) becomeswithBy the Routh–Hurwitz Criterion, all roots of the polynomial (26) have negative real parts if and only if , , and . When , we have and . Since , we only need to consider .
After some algebraic manipulations ([21]) we can prove thatSo we conclude that when and , the endemic equilibrium is locally asymptotically stable for .
Consider now the case when 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 be a purely imaginary roots of (24), thenwhich is equivalent toSeparating real and imaginary parts leads toAdding the squares of both equations together givesUsing Lemma 1 and after some algebraic manipulations, equation (32) can also be written in the following form:where are as in Theorem 2. If , 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 , which can be considered as a perspective of this work. Consider , then (24) becomesif is a root of (34). After substituting and separating real and imaginary parts, we have
Adding the squares of both equations together giveswhere , , and are as follows:
Denote the biggest positive root of (37); then from (35), we have
Let
Then, we can define as the first value of when characteristic roots cross the imaginary axis.
Further, differentiating equation (34) with respect to , we get
This givesand after some algebraic manipulations, we get
Thus,and
So, if and , there exists a Hopf bifurcation as .
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 , 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 less than 1, and we shall do the same for parameters verifying 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 corresponds to cells [22]. For our numerical simulation, we consider cell populations x and y, virus population v is expressed in units of , 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 .
Figure 2 presents the curves of system (2), using various initial conditions when and . In Section 3, using a suitable Lyapunov function, we have proved that, in this case, , the disease-free equilibrium is globally asymptotically stable. From this figure, we see that the curves converge to the free equilibrium , 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 , , 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 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 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 is asymptotically stable and, for , the equilibrium point is unstable, and when , a Hopf bifurcation of periodic solutions of system (2) occurs at (Figure 7).
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 ratio , 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 , and when , the local stability of the endemic equilibrium point depends on function , 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 , the virotherapy fails as the population of tumors cells increases and the population of infected tumor decreases, and if , 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].
Data Availability
Te disciplinary data used to support the findings of this study have been deposited in the Network Repository (http://www.networkrepository.com).
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.