Backward Hopf bifurcation in a mathematical model for oncolytic virotherapy with the infection delay and innate immune effects

In this paper, we consider a system of delay differential equations that models the oncolytic virotherapy on solid tumours with the delay of viral infection in the presence of the innate immune response. We conduct qualitative and numerical analysis, and provide possible medical implications for our results. The system has four equilibrium solutions. Fixed point analysis indicates that increasing the burst size and infection rate of the viruses has positive contribution to the therapy. However, increasing the immune killing infection rate, the immune stimulation rate, or the immune killing virus rate may lead the treatment failed. The viral infection time delay induces backward Hopf bifurcations, which means that the therapy may fail before time delay increases passing through a Hopf bifurcation. The parameter analysis also shows how saddle-node and Hopf bifurcations occur as viral burst size and other parameters vary, which yields further insights into the dynamics of the virotherapy.


Introduction
Oncolytic virotherapy is an encouraging therapeutic strategy to destruct tumours [6]. Oncolytic viruses selectively infect and replicate in tumour cells, but spare normal cells. It was initially tested in the middle of the last century, and merged with renewed interests over last 30 years due to the technological advances in virology and in the use of viruses as vectors for gene transfer [14]. However, the efficacy of oncolytic viral therapy has not been well established yet. The major challenge is the complexity of the immune responses [5,38] in the process of the therapy [7].
Since fifteen years ago, mathematicians have applied mathematical models to understanding of oncolytic virotherapy. Wu [37] and Wein [33] proposed and analysed a system of partial differential equations (PDEs) that is essentially a radially-symmetric epidemic model embedded in a Stefan problem to describe some aspect of cancer virotherapy.
Wodarz [34,35] formulated a simple model with three ordinary differential equations (ODEs) to study three hypothetical situations, viral cytotoxicity kills tumour cells, a virusspecific lytic CTL response contributes to killing of infected tumour cells, and the virus elicits immunostimulatory signals within the tumour that promote the development of tumour-specific CTL. Komarova [16] and Wodarz et al. [36] analysed several possible mathematical formulations of oncolytic virus infection in terms of two ODEs, while Novozhilov et al. [19] analysed ratio based oncolytic virus infection models. Bajzer et al. [3] used three ODEs to model specific cancer virotherapy with measles virus, and then considered optimization of viral doses, number of doses, and timing [4].
Friedman et al. [12] proposed a free boundary problem to study brain tumour glioma with mutant herpes simplex virus therapy. It incorporated immunosuppressive agent cyclophosphamide to reduce the effect of the innate immune response. This model reveals the oscillation phenomenon in oncolytic viral therapy. Vasiliu et al. [30] studied a simplified model for the cell population oscillation, which may be caused by interaction between infected tumour cells and innate immune cells. To obtain a basic dynamic picture of oncolytic viral therapy, Tian [27] proposed a simple model with three ODEs to represent interaction among tumour cells, infected tumour cells, and oncolytic viruses, and concluded that the viral therapeutic dynamics is largely determined by the viral burst size. To further understand how the viral burst size is affected, Wang et al. [31] and Tian et al. [28] proposed a delay parameter for virus infection process into the basic model. These delay differential equation models gave another explanation of cell population oscillation and revealed a functional relation between the virus burst size and lytic cycle. In a recent paper [8], Choudhury et al. considered a simple model of three ODEs for the dynamics of oncolytic virotherapy in the presence of immune response. However, this model did not include the free virus population, and it may not give a complete the dynamic picture of viral therapy with innate immune response. There also are some other mathematical models that describe the dynamics of oncolytic viral therapy [20,24].
All proposed mathematical models have given some insights into oncolytic virotherapy. However, there is a considerable need of understanding the dynamics of oncolytic virotherapy in the presence of immune responses [7]. Particularly, we need to understand the subtle dynamics which includes some details of viral infection process and immune responses. The viral oncolytic process has five stages, attachment, penetration, replication, assembly, and release [18], and each stage is preprogrammed and tightly regulated in time and space [25]. There are several articles which give a detailed description about the five stages of viral infection, for example, the reader is referred to review article [21], [24] and [39] for more details. To infect a tumour cell, one or more virus particles must enter the host cell. This process differs between classes of viruses, especially regarding the presence or absence of a viral membrane envelope. Attachment to cells is mediated by viral surface proteins targeting molecules accessible on cell membranes. Upon attachment, the viral particle passes the host cell membrane in a process termed penetration, which is involved in different mechanisms depending on the virus type. Some viruses require interaction of two distinct surface proteins for entry [21]. Binding of the H protein to a cell entry receptor induces a change in the structural conformation of the F (fusion) protein, leading to an approximation of viral particle and host cell. This finally results in membrane fusion, allowing entry of the viral particle into the cell [17]. Endocytosis provides an alternative route of cell entry. Acidification of maturing endosomes can induce conformation changes in viral attachment proteins to promote membrane approximation and fusion [26]. Nonenveloped viruses do not need for membrane fusion and rely on endocytic pathways for entry [39]. After one or more viruses enter a host cell, the virus will replicate up to certain number of copies, then assemble, and then release from the host cell by lysis. It is clear that each stage needs some time to complete [18]. In order to understand the subtle dynamical behaviours of the oncolytic viral therapy, it is necessary to incorporate each stage of viral oncolytic process into mathematical modelling. One way to include these infection stages into a mathematical model is to introduce delay parameters for each stage. However, it is reasonable to use two parameters to present the viral oncolytic process, the infection time which is the time period from the time when viruses attach to a host cell to the time when they penetrate into the host cell, and the viral renewal time (or the lytic cycle in narrow sense) which is the time period from the time when viruses penetrate to the host cell to the time when new viruses burst out from the host cell. In the literature, there are several models which include the infection time as a delay parameter [28,31,32], and several models which include the virus renewal time as a delay parameter [1,2,15]. However, these models either have two equations or three equations for tumour cells and infected tumour cells or/and viruses, and do not incorporate immune responses. To obtain insights, we have built a mathematical model for oncolytic virotherapy with the innate immune system [20], and a computational model [29]. In the current study, we further build the infection time and innate immune response into our basic model as the system of delay differential equations that models the oncolytic virotherapy on solid tumours with the delay of viral infection in the presence of the innate immune response.
Our mathematical model incorporates tumour cells, free viruses, infected tumour cells, and innate immune cells. The infection time of oncolytic viruses is considered as a time delay parameter, which varies from seconds to hours [9,18,22]. We conduct analysis and numerical computations. Our model shows some interesting dynamical properties. For example, our model has saddle node bifurcations and Hopf bifurcations as different parameters vary. Biologically, our analysis shows that if we increase the burst size of the oncolytic virus and infection rate, we can obtain good therapeutic results. If we increase the immune killing infection rate, the immune stimulation rate, or the immune killing virus rate, the oncolytic viral therapy may fail. The infection time delay induces backward Hopf bifurcations, which means that we will have unstable periodic solutions before time delay increases passing a Hopf bifurcation.
The rest of the paper is organized as follows. In Section 2, we provide basic feasibility of the model. In Section 3, we analyse and simulate our model without infection delay. In Section 4, we analyse and simulate our model with infection delay. We then conclude our study with Summary and discussion in Section 5.

Mathematical model and positiveness
We assume the solid tumour follows an exponential growth pattern as we are concerned with glioma viral therapy [12], and as in the early stage of tumour initiation the tumour is in a fast growth phase. We only consider the infection time as a delay parameter. Let x(t), y(t), v(t), and z(t) represent the number of tumour cells, infected tumour cells, free virus and innate immune cells, respectively. Based on our models [20,31], the mathematical model for oncolytic viral therapy with the infection delay and the innate immune response is proposed in the following. (1) The parameters are defined as follows. The tumour growth rate is a. β is the viral infection rate. τ is the infection time delay. μ is the immune killing infection rate that is the rate at which the innate immune system kills the infected tumour cells. δ is the lysis rate of infected cells. b is a key parameter standing for the virus burst size. k is the immune killing virus rate that is the rate at which the innate immune system kills free viruses. γ is the clearance rate of viruses. r is the immune stimulation rate that is the rate at which the infected tumour cells stimulate innate immune response. c is the constant determining the clearance of immune cells. The infection, lysis, and the mechanism of virotherapy are shown in Figure 1.

Model analysis in the absence of infection delay
In this section, we study the existence and stability of equilibria of the model (1). We numerically analyse how a saddle node bifurcation occurs as main parameters change. We also numerically analyse how a Hopf bifurcation occurs when the burst size b varies.

Analytical results
The system (1) has two equilibria on the boundary as follows.
components of E 1 are positive, and it is unstable since there is at least one positive eigenvalue.
For the possible existence of positive equilibria of model (1), we calculate and find there are at most two inner equilibria given by E ± : (x ± , y ± , v ± , z ± ) with v ± = a/β, y ± = (c/r)z ± , x ± = (bδc/ar − k/β)z ± − γ /β and z ± are the roots of For definiteness, we always assume z + ≥ z − . Since x ± is monotonically increasing in z ± , then we have x + ≥ x − . Obviously b = b * is a saddle-node singularity. We now mainly focus on the stability of E ± .
When τ = 0, the Jacobean matrix of the right hand side of the system (1) is given by At each of the equilibria E ± , we can compute this matrix by To determine the stability of E ± , we need the characteristic equation of J(x ± , y ± , v ± , z ± ) which is where, after dropping all subscripts to simplify notations, βr .
The stability of E + can be easily determined. In fact, we have

Lemma 3.2:
which is negative at z = z − from (2). If z = z + , then obviously we have A 42 > 0, which means A 4 < 0. This indicates the product of all four roots of the characteristic equation is negative, which means there exists at least one positive root, i.e. E + is always unstable.
Noticing that 1 = A 1 > 0 and using the Routh-Hurwitz criterion, we have These two Lemmas are conditions to determine the stability of E ± , and also indicate that only Hopf bifurcation can lead to instability of E − . In fact, we may vary some parameter values, for example, increase b, so that i.e. the Routh-Hurwitz condition is violated, then a possible Hopf bifurcation occurs at b = b H .

Numerical simulation
It seems difficult to obtain analytical results about the occurrence of Hopf bifurcations. However, we can numerically demonstrate this. Choose a set of parameter values as in [12], Using these parameters we can compute that E − = (26105.7, 994.0, 100000, 9.9). The four key values determining the stability of E − are, respectively, They are all positive, hence Theorem 3.1 yields that E − is stable when τ = 0, which is shown in Figure 2.
Let b, μ, r, k and β vary separately, we have the saddle-node bifurcation diagram shown in Figure 3. From these computations, we find that, if we increase the burst size b or infection rate β, the system moves to a stable E − . This means that the viruses control the tumour growth and the therapy succeeds. If we increase immune killing infection rate μ, immune killing virus rate k, or immune stimulation rate r, the viruses cannot control the tumour growth and the therapy fails.   Let the parameter b vary, and keep other parameters as in (P). We can calculate the values of 3 as a function of b, and its graph is shown in Figure 4. In this figure, we find when b H = 153.649, 3 = 0. We also see that, when parameters fall in the left of the graph, we have stable equilibria; when parameters fall in the right of the graph, we have periodic solutions induced by Hopf bifurcations. The numerical simulations in Figures 5 and 6 indicates there exist stable Hopf bifurcating periodic solutions.

The effect of delay
In this section, we study the effect of delay τ at the equilibrium E − , and Hopf bifurcation induced by this infection delay. Since the system has four state variables and the research about the characteristic equation is difficult to be carried out analytically, we will combine numerical analysis in this section to show the dynamical behaviour of the system with infection delay.

Analytical results
In the presence of time delay, the characteristic equation at E − is given by where all the parameters B j and C j have quite complicated expressions, listed as follows, and When τ = 0, we assume all roots of (3) have negative real parts. That is, the Routh-Hurwitz condition in Theorem 3.1 holds true. As discussed in the previous section, we only consider the case of Hopf bifurcations.
For τ > 0, plugging λ = iω into the characteristic equation (3) yields which gives Obviously, at most four positive roots of ω can be solved from (5), which are denoted by ω j , j = 1, 2, 3, 4. To ensure the existence of Hopf bifurcation we assume the transversality condition holds true, that is Re(dλ/dτ ) = 0 at τ = τ j l for any l and j. Now we have the main theorem about delay inducing Hopf bifurcation [23]. When the parameter value τ increases, Hopf bifurcation occurs at τ 0 , τ 1 , . . .. E − is locally stable for all τ < τ 0 . Thus the first Hopf bifurcation value τ 0 is important for the viral therapy.
In fact, we can calculate the stability of the bifurcating periodic solutions and the direction of the Hopf bifurcation on the local centre manifold [13]. The process is standard but with very complicated and tedious, thus we only give some numerical results in the next subsection instead of proceeding analytically.
When τ < τ 0 , E − is locally stable. We take τ = 80 which is less than τ 0 , as shown in Figure 7, the locally stable E − and an unstable periodic solution coexist. Figures 8 and 9 show that solutions with small (or large) initial values converges to E − (diverges to infinity).
When τ > τ 0 , from Theorem 4.1, the equilibrium E − is locally unstable, and the figure is not shown here.
It is better to study how parameters are related to each other for the bifurcations. In Figure 10, we show the bifurcation diagrams on parametric planes τ 0 − b, τ 0 − δ, τ 0 − r, τ 0 − γ , τ 0 − μ, and b − δ. In these figures, we calculate the Hopf bifurcation values.   E − loses its stability when τ > τ 0 and the Hopf bifurcation occurs. The oncolytic virotherapy needs to maintain the stability of E − in order to succeed. From the figures in 10, we can see how we can avoid bifurcation occurrence in terms of parameter values. For the lysis size b and τ 0 , if τ 0 is small, we may take b big. For the infected tumour cell burst rate δ and τ 0 , we have a similar qualitative conclusion as for the burst size b and τ 0 . For the immune killing infection rate μ and τ 0 , if τ 0 is small, we have to take μ small. For the immune stimulation rate r and τ 0 , or the virus clearance rate and τ 0 , we have a similar qualitative conclusion as for the immune stimulation rate r and τ 0 . One interesting observation is that, when the infection time τ is fixed, if the infected tumour cell lysis rate δ is small, then the burst size can be big.

Summary and discussion
In this research, we study a 4-dimensional delay differential equation system, which models oncolytic viral therapy with the innate immune response and infection delay. We analyse and numerically simulate our model. We obtain that, if we increase the viral burst size and infection rate of the virus, it will have positive contributions to the therapy; if we increase the immune killing infection rate, the immune stimulation rate, or the immune killing virus rate, the therapy may fail. These conclusions seem reasonable. The viral burst size is a measure of virus replication. If the viral burst is big, the tumour site will have more viruses which increase the chance of infection. The virus infection rate is also a measure of ability for which the virus can infect tumour cells. If the infection rate is big, the tumour cells will have high probability to be infected. These both will increase the probability of success of the therapy. However, if we increase the immune killing infection rate, the immune stimulation rate, or the immune killing virus rate, we decrease the infection (infected tumour cells), get more innate immune cells that decreases infected tumour cells potentially, or decrease the number of free viruses. Overall, increasing these three parameters decrease the infection of the tumour.
The viral infection time delay induces backward Hopf bifurcations, which is different from the forward Hopf bifurcation induced by the viral infection time delay in virotherapy model without the innate immune response [31]. On the other hand, the backward Hopf bifurcation means that there will be unstable periodic solutions before time delay passes through a Hopf bifurcation. We also find that the Hopf bifurcation value increases while the burst size decreases. This may suggest that when biologists change virus genome for a certain burst size, we should consider its effect on virus infection time if we want a stable outcome of the therapy. Since the Hopf bifurcation induced by the infection time delay is backward, it has an effect that the therapy may fail before time delay increases passing through a Hopf bifurcation.
It is clear that we could eradicate tumours if there would be only such simple interactions with the immune system. The fact is that there are not only the innate immune response but also the adaptive immune response in the oncolytic viral therapy. And more complicated, these two immune responses may play opposite roles in viral therapies. Furthermore, there are several steps in the viral infection process, particularly, the infection time and virus renewal time, which may substantially influence the whole oncolytic outcome. This calls for further mathematical modelling to find the balance between the two immune responses with the infection process. We are now working on this direction in our next models.