Mathematical model of cancer treatment with virotherapy and immune

system

Arvind Kumar Sinha1 , Ambika1∗

National Institute of Technology Raipur, Department of Mathematics, Raipur, Chhattisgarh, India

1
Arvind Kumar Sinha, aksinha.maths@nitrr.ac.in, ORCID ID: 0000-0002-8684-724X
1∗
Corresponding Author: Ambika, ambika11sahu@gmail.com

Abstract. Cancer affects people worldwide. Although traditional treatments like chemotherapy and radiation are
effective, they have side effects and harm healthy cells. Virotherapy is a treatment that specifically targets cancer
cells, leaving healthy cells unharmed. This paper presents a mathematical model to explore the dynamic interplay
among cancer, virotherapy, and the immune response. The model considers virus replication within the tumor, im-
mune cell infiltration, and the tumor-immune microenvironment. We verify the model by showing that the solutions
are non-negative and bounded. Our model contains five equilibrium points: dead equilibrium point, cancer invasion
equilibrium point, immune response-free equilibrium point, virus-free equilibrium point, and coexisting equilib-
rium point. We find the local stability and global stability of the equilibrium points. We simulate numerically to
demonstrate our analytical results. A high stimulation immune response rate and a large burst size reduce tumor
size. The paper presents the importance of high immune system and virotherapy on cancer patients. The model
helps to oncologist and immunologists to find the behavior of virotherapy and immune response to control the
proliferation of different kinds of tumors.

Keywords: Mathematical model, Tumor, Virotherapy, Immune response, Stability

1 Introduction

Cancer tumors are dangerous for humankind, and they may cause death. In 2020, nearly 10 million deaths by cancer
caused, according to World Health Organization (WHO)[1]. Efforts are underway to discover new cancer treatments
due to the low efficacy and high toxicity of traditional therapies like surgery, chemotherapy, and radiation. Therefore,
efforts to find new cancer treatments are a nonstop process. Virotherapy is a new treatment made possible by advances
in genetic engineering that researchers have extensively studied. Virotherapy is a treatment that uses genetically modi-
fied viruses to infect and replicate within cancer cells (CC) while leaving normal cells unharmed. Virotherapy appears
to be a favorable treatment but faces the challenge of evading immune cells (IC) that destroy viruses before they infect
CC [2]. Clinical trials have demonstrated the effectiveness of oncolytic viruses (OV) [3].
An immune response (IR) is a physiological reaction that occurs within an organism to defend against external factors
in the context of inflammation [4]. However, the immune response plays a vital role in cancer cells. More recently, it
has become apparent that innate immune responses also play an essential role in resistance against tumor development
and progressive growth. The immune response presents a challenge in maximizing efficacy.
Various mathematical models have been developed to understand the interaction between CC and OV regarding the
dynamics of virotherapy. Wodarz [5] proposed a fundamental model to describe the growth of tumors when treated
with virotherapy. The model shows how uninfected cancer cells (UCC), infected cancer cells (ICC), and free viruses
interact. Bajzer et al. [6] proposed a model that eradicates free viruses from cancer cell infection. They also introduced
a generalized logistic growth model for cancer cells. Their model provides an explanation for the direct transmission
of infection from infected cells to uninfected cells. Tian’s model in [7] shows that the virus burst size is critical for in-
fection spread. Tian’s model indicates that larger burst sizes can decrease tumor cells (TC). Wodarz [8] conducted one
of the initial studies on the impact of IR on tumor virus dynamics. He proposed a model that describes the dynamics
among the tumor, the virus, and virus-specific Cytotoxic T Lymphocytes (CTL). The model comprises three variables:
uninfected tumor cells, tumor cells infected by the virus, and virus specific CTL. Ashyani et al. [9] analyzed a modi-
fied model [8] with simulations that combined IR to virus and tumor as a single variable. Dependent variables of the
models in [8,9] are UCC, ICC, and IC. Phan and Tian [10] introduced an additional state variable to describe the vi-
rotherapy. They analyzed the impact of the innate IR on ICC and virotherapy. After analyzing the data, they concluded
that the model behaves similarly to [7] without innate immunity when the burst size is large. When the burst size is
small, innate immunity produces more equilibria in the dynamics. Tuwairqi et al. [11] proposed a modified model
based on Phan and Tian’s [10] model. The modified model includes UCC, ICC, virotherapy, and an immune response.
They have added three new parameters to the previous model. In prior models, the stability of virotherapy outcomes
was investigated through analytical and numerical methods, presenting different criteria for partial reduction, tumor
eradication, and therapy failure. Furthermore, these models only consider the variations in cell densities concerning
time. The papers [12,13] present research on the spatio-temporal dynamics of TC that undergo virotherapy and radio
virotherapy treatments. These studies involve the use of partial differential equations (PDEs). General biological mod-
els involving cell density diffusion are presented in [14,15,16].
In this paper, we introduce a more realistic interaction model among tumor cells, virotherapy, and immune system. It
is based on the fact that the immune system functions as a predator, which works like a natural killer cell. Our model
incorporates immune response stimulation rates and tumor destruction via cytokine release from natural killer cells.
Thus, the dynamical analysis identifies the optimal treatment and appropriate options for achieving cancer remission.
The whole work is divided into 6 sections. Section 1 is the introduction part of the work then the next section 2 presents
the mathematical model of the proposed dynamics. In section 3, we have analyzed the model in three sub-sections.
Section 3.1 deals with positivity and boundedness, section 3.2 calculates the critical points, and section 3.3 studies
the stability of the critiacl points. In section 4, we illustrate the numerical parameter analysis of model with analytical
results and numerical experiments. Section 5 includes the model discussion, while section 6 includes the conclusions.

2 Mathematical model

This model explains the interaction among TC, virotherapy, and immune response (IR). It is composed of three groups:
cancer cells T (t) , virotherapy V (t) , and immune response I(t). CC grows logistically with r growth rate and k is
carrying capacity and dies at a rate of d. OV infects CC at a rate of a. IC kill CC at a rate of b. e is burst size . IC
eliminate the virus with a rate of c. Virus decay at a rate f . CC activate the IR at a rate m. Eventually, the IC decay at
a rate n.
The model is governed by the following systems of non-linear equations :

2
 
dT

 = rT (1 − kT ) − aT V − bT I − dT,
dt







dV

= eT − aT V − cV I − f V, (2.1)
 dt





 dI = mT I − nI,



dt
with initial conditions

T (0) = T0 ≥ 0, V (0) = V0 ≥ 0, I(0) = I0 ≥ 0.

3 Qualitative analysis

We examine the model in this section through qualitative analysis. First, we define the positivity and boundedness
regions for the dependent variables. Then, we get and analyze the stability of model constant solutions.

3.1 Non-negativity and boundedness

n 1
Theorem 1. Solutions of the model are non-negative and bounded in region Γ1 = (T, V, I) ∈ R3+ | T ≤ , V ≤
k
e b r o
,T + I ≤ , with non-negative initial conditions in the same region Γ1 .
f m σ1
Proof. From the model
dT
= ϕ1 (T, V, I) dt, where
T
ϕ1 (T, V, I) = r(1 − kT ) − aV − bI − d. Integrating over [0, t], gives
hZ t i
T (t) = T (0) exp ϕ1 (T, V, I)ds
0

∵ T (0) ≥ 0, we have T (t) ≥ 0 for all t ≥ 0.

Also, from the model

dV
≥ −(aT + cI + f )V, then it can be written as
dt
dV
≥ ϕ2 (T, I)dt, where ϕ2 (T, I) = −(aT + cI + f ). Integration over [0, t] gives
V
hZ t i
V (t) ≥ V (0) exp ϕ2 (T, I)ds .If V (0) is greater than or equal to 0, then we have V (t) ≥ 0 for all t ≥ 0.
0

And also again from the model

dI
= ϕ3 (T )dt, where ϕ3 (T ) = (mT − n).
I
Therefore,
hZ t i
I(t) = I(0) exp ϕ3 (T )ds . Also, since I(0) is greater than or equal to 0, we have I(t) ≥ 0 for all t ≥ 0. Hence,
0

3
(T (t), V (t), I(t)) is non-negative.

For boundedness: we begin from the model

dT
≤ rT (1 − kT ).
dt
du
Now, we assume that the differential equation = ru(1 − ku) with initial condition u(0) = u0 . We check for the
dt
following solution:
u0 1
u(t) = . Hence, lim sup u(t) = .
u0 k + (1 − u0 k)e−rt t→∞ k
dT dV 1
We know ≤ . Thus, lim sup T (t) ≤ lim sup u(t). Therefore, lim sup T (t) ≤ .
dt dt t→∞ t→∞ t→∞ k
Again, from the model

dV
≤ eT − f V
dt
e
≤ e − f V, we obtain lim sup T (t) ≤ . Let
t→∞ f
b dW dT b dI
W (t) = T (t) + I(t). Then = +
m dt dt m dt
b b
= rT (1 − kT ) − bT I − dT + mT I − nI
m m
b
≤ rT − dT − nI
m
b
≤ r − σ1 (T + I),
m
where σ1 = min(d, n).
r
Thus, lim sup W (t) ≤ .
t→∞ σ1

Hence, solutions of the model are non-negative and bounded in Γ1 :

n 1 e b r o
Γ1 = (T, V, I) ∈ R3+ | T ≤ , V ≤ , T + I ≤ .
k f m σ1
Γ1 is positively invariant, which means that any solution starting in Γ1 remains there ∀ t ≥ 0.

3.2 Equilibrium (Critical) points

The system’s Critical points are constant solutions found by setting the rates in equation (2.1) to zero. We have found
five critical points Ei∗ = (Ti , Vi , Ii ), where i = 0, 1, 2, 3, 4 :
1 d   (r − d)(a + f ) e(r − d)  n rm − (rkn + dm) 
E0∗ = (0, 0, 0), E1∗ = 1− , 0, 0 , E2∗ = , , 0 , E3∗ = , 0, ,
k r rk(a + f ) + ae rk(a + f ) + ae m bm
E4∗ = (T4 , V4 , I4 ),
where E4∗ satisfies the following equations:
r(1 − kT4 ) − aV4 − bI4 − d = 0,
eT4
V4 = ,
aT4 + cI4 + f
mT4 − n = 0.
E0∗ denotes dead critical point where all tumor cells, viruses and immune cells tend to zero. E1∗ denotes cancer invasion

4
equilibrium point, E2∗ denotes immune response-free equilibrium point, E3∗ denotes a fourth equilibrium point where
viruses tend to zero. As for E4∗ , denotes coexisting critical point.

3.3 Stability of equilibrium

We use the linearization method [17] to assess the local stability (LS) of the critical points. The method involves eval-
uating the Jacobian matrix at the equilibrium point and proving that the eigen values of the linearized system’s matrix
are negative. Moreover, we use Lyapunov functions to examine the global stability (GS) of critical points [18].

For locally asymptotically stable (LAS)

Theorem 2. If r < d , then E0∗ is LAS.

Proof. The Jacobian matrix of system at E0∗ :

 
r−d 0 0
 
J(E0∗ ) = 
 e −f 0 

0 0 −n

Thus, the eigen values of J(E0∗ ) are λ1 = r − d, λ2 = −f, λ3 = −n.

All eigen values are negative, and λ1 is negative if r < d. Hence, E0∗ is LAS if r < d.

m(r − d)
Theorem 3. If < 1 , then E1∗ is LAS.
nkr
Proof. System’s Jacobian matrix at E1∗ , we get

−a d −b d
 
d−r (1 − ) (1 − )
 k r k r 
∗
 −a d 
J(E1 ) =  e
 (1 − ) − f 0 
k r 
 m d 
0 0 (1 − ) − n
k r
We solve the characteristic equation |J(E1∗ ) − λI|= 0, then we get the eigen values :
m d
λ1 = (1 − ) − n and the equation λ2 + a1 λ + a2 = 0,
k r
a d ar d ae d ad d
where a1 = (1 − + f + r − d), a2 = (1 − ) + rf + (1 − ) − f d − (1 − ).
k r k r k r k r
m(r − d)
It is obvious that λ1 is negative if < 1. By using Routh Hurwitz criteria (RHC), the eigen values
nkr
λ2 and λ3 are negative if a1 > 0 and a2 > 0. Clearly a1 and a2 are positive. Hence, E1∗ is LAS if
m(r − d)
< 1.
nkr
Theorem 4. E2∗ is LAS if mT2 < n and e > aV2 .

5
Proof. System’s Jacobian matrix at E2∗ is

 
−rkT2 −aT2 −bT2
 
J(E2∗ ) = 
e − aV2 −aT2 − f −cV2  
0 0 mT2 − n

After solving the characteristic equation |J(E2∗ ) − λI|= 0, we find one eigen value explicity, λ1 = mT2 − n,
Whereas other eigenvalues λ2 ,3 satisfy the following equation: λ2 + a1 λ + a0 = 0,
where
a1 = rkT2 + aT2 + f, a0 = arkT22 + rf kT2 + aeT2 − a2 T2 V2 .

It is clear that a1 > 0. As for a0 , we have

a0 = arkT22 + rf kT2 + aeT2 − a2 T2 V2

= arkT22 + rf kT2 + aT2 (e − aV2 ).

Thus, for a0 to be positive we must have e > aV2 . Hence, E2∗ is LAS if mT2 < n and e > aV2 .

Theorem 5. E3∗ is LAS provided ai > 0 where i = 0, 1, 2, and a2 a1 > a0 . Here a0 , a1 and a2 are provided in the
proof.

Proof. System’s Jacobian matrix at E3∗ , we get

 
−rkT3 −aT3 −bT3
 
J(E3∗ ) = 
 e −aT3 − cI3 − f 0 


mI3 0 0

The characteristic equation |J(E3∗ ) − λI|= 0 gives λ3 + a2 λ2 + a1 λ + a0 ,

where
a2 = rkT3 +aT3 +cI3 +f, a1 = rakT32 +rkT3 z3 +rkT3 f +aT3 e+bmT3 I3 , a0 = bT32 amI3 +bmcT3 I32 +bT3 mI3 f .

It is clear that a2 > 0, a1 > 0 and a0 > 0. We need to show that a2 a1 > a0 , then by using RHC, we prove that
eigen values λ2 ,3 are negative.
Now,
a2 a1 − a0 = ar2 k 2 T33 + r2 k 2 T32 I3 + r2 k 2 f T32 + arkeT32 + brmkT32 I3 + a2 rkT33 + arkT32 I3 + arkf T32 + a2 eT32 +
arckT32 I3 + rckT3 I32 + rckf T3 I3 + aceT3 I3 + arkf T32 + rkf T3 I3 + rkf 2 T3 + af eT3 .

Thus, E3∗ is LAS if a2 a1 − a0 > 0.

Theorem 6. E4∗ is LAS provided ai > 0, where i = 0, 1, 2, and a2 a1 > a0 . If f > (r + mT4 ).

6
Proof. Evaluating the Jacobian at the critical point E4∗ and solving the characteristic equation |J(E4∗ ) − λI|= 0, we
get the equation λ3 + a2 λ2 + a1 λ + a0 = 0,
where
a2 = 2rT4 k − (r + mT4 ) + aV4 + bI4 + d + aT4 + cI4 + f + n

a1 = mT4 r − r(aT4 + cI4 + f + n) + 2rT4 k(aT4 + cI4 + f + n) + aV4 (aT4 + cI4 + f + n) + bI4 (aT4 +
cz4 +f +n+mT4 )−mT4 (aT4 +cI4 +f )−mT4 (bI4 +d+2rT4 k)+d(aT4 +cI4 +f +n)+ncI4 +f n−a2 T4 V2 +aT4 e

a0 = mT4 r(aT4 + cI4 + f ) − nr(aT4 + cI4 + f ) + 2rT4 k(nT4 + f n) − 2rT 2 k(mT4 + cI4 + f m) + a2 V4 T4 (n −
mT4 ) + aV4 I4 (ncm T4 c) + aV4 f (n − mT4 ) + bI4 n(aT4 + cI4 + f ) − bI4 mT4 (aT4 + cI4 + f ) − mT4 d(aT4 + cI4 +
f ) + nd(aT4 + cI4 + f ) − aT4 (meT4 + aV4 n + mcV4 I4 ) + aenT4 + a2 mT 2 V + bmT4 I4 (aT4 + cI4 + f ).

Clearly, a2 > 0 if f > (r + mT4 ) . Hence, by the RHC, the LAS of E4∗ which is conditions stated in the theorem.

For Global asymptotically stability (GAS)

Theorem 7. E0∗ is GAS if r < d.

Proof. Define the Lyapunov function as

L1 (T, V, I) = mcT + V + bcI.

It is evident that L1 is a positive definite function. By computing the derivative of L1 along the solutions of the model,
we obtain

L′1 = mc(rT (1 − kT ) − aT V − bT I − dT ) + eT − aT V − cV I − f y + bc(mT I − nI).

= mc(r − d)T − mcrT 2 k + T I(−mbc + mbc) + eT − aT V − cV I − f V − bcnI
≤ mc(r − d)T − f V − bcnI.

If r < d, then L′1 ≤ 0 and L′1 = 0 iff T = V = I = 0. Thus, the largest invariant set in (T, V, I) ∈ Γ1 : L′1 = 0 is E0 .
Hence, by La Salle’s invariance principle, E0∗ is GAS.

fr 1 cmrI12 2
Theorem 8. E1∗ is GS if 1 − < 0 and R0 < 1, where R0 = (T1 + ) .
a rI1 cn 4
Proof. Define the Lyapunov function,

2
L1 (T, V, I) = (T − T1 )2 + V + bcI.
rT12

7
 By differentiating L1 and using d = r(1 − kT1 ), we get

4
L′1 = (T − T1 )T ′ + V ′ + bcI ′
rT12
4
= (T − T1 )(−krT 2 − aT V − bT I + krT T1 ) + eT − aT V − cV I − f V + bc(mT I − nI)
rT12
4 T1 T2
= 2 (−krT (T − T1 )2 − aV (T − )2 + aV 1 − bT 2 I + bT T1 z) + eT − cV I − f V + bc(mT I − nI)
rT1 2 4
Collecting and simplifying terms, we obtain

a −4 4
L′1 ≤ V ( − f ) + I( 2 bT 2 + ( bT1 + bcm − bcn)T )
r rT1 rT1
a 4 1 bcmrT12 2 1 bcmrT12 2
= V ( − f) − bI(T − (T 1 + )) + bI(T1 + ) − bcnI
r rT12 2 4b rT12 4b
a 1 bcmrT12 2
≤ V ( − f) + bI(T1 + ) − bcnI
r rT12 4b

fr
If 1 − < 0 and R0 < 1, then E1∗ is GS.
a

4 Numerical analysis

We numerically solve the system in this section and present various simulations of the model using the fourth-order
Runge-Kutta method. We analyze parameter sensitivity to identify critical ones for optimal therapeutic strategies. The
values of the parameters are taken from[10,19,20,21] and are shown in Table 1.

Table 1. Parameters description

Parameter Description Value Unit Reference

−2
r logistic growth rate of CC 2 × 10 1/h [19]
k Carrying capacity of the environment for CC 1.5 1/h estimated
d Natural death rate of CC 0.0071 1/h [20]
−9 3
a The rate at which CC infected 7/10 × 10 mm /h virus [19]
and killed by the OV
b Immune killing rate of CC 2 × 10−8 mm3 /h immune cell [11]
e Burst size of free virus 50 virus/cell [19]
−8 3
c Immune killing rate of the virus 10 mm /h immune cell [19]
f Clearance rate the virus 0.0119 1/h [21]
−7 3
m Stimulation rate of the IR by CC 5.6 × 10 mm /h cancer cell [11]
−3
n Clearance rate of the IR 2 × 10 1/h [10]

8
4.1 Numerical experiments and Parameter analysis

We solve the system numerically with initial values:

T (0) = 0.5, V (0) = 0.01, I(0) = 0.01.
We change some parameter values in Table 1 to ensure stability for each critical point. First, we change the value
of k = 1.5, d = 0.3, e = 0.05, f = .99, and n = 1. Figure 1 shows that E0∗ = (0, 0, 0) is LAS, if r < d. If
we change the parameters r = 0.4, k = 1.4, d = 0.1, e = 0.001, f = 0.99, and n = 2. Figure 2 shows that
m(r − d)
E1∗ = (0.5357, 0, 0) is LAS, if < 1. If we change the parameter r = 0.45, k = 2, d = 0.23, e = 0.999, f =
nkr
0.19, m = 0.0005, and n = 1, the conditions for the presence and stability of E2∗ are met. Therefore, figure 3
shows the LAS of E2∗ = (0.2529, 1.3084, 0). E3∗ = (0.0916, 0, 0.8278) is LAS as seen in figure 4. If we change
the parameters: r = 0.4, k = 2, d = 0.1, e = 0.315, f = 0.69, and n = 0.0001, then figure 5 shows the LAS of
E4∗ = (0.3759, 0.1791, 0.1000).
Parameter analysis: Here, we explore changes in the model’s dynamics caused by varying burst size (e) and the
immune response stimulation rate by cancer cells (m). We establish the initial values as follows: (T = 0.5, V =
0.01, and I = 0.01).
The burst size (e) defines the number of new viruses released from the lysis of an infected CC. Fig. 6 displays how
cancer cells change with different burst sizes over time.
As for the relative immune response stimulation rate, m. Fig. 7 shows the time variation on the cancer cell for different
values of stimulation immune response.

Fig. 1. Time variation in the relative population sizes with initial values. Here parameter values are: r = 0.02, a = 7/10×10−9 , b =
2 × 10−8 , k = 1.5, d = 0.3, e = 0.05, c = 10−8 , m = 5.6 × 10−9 , and n = 1

5 Discussion

Virotherapy is one of the type of cancer treatment that uses genetically modified viruses to target and destroy CC
and is a promising cancer treatment approach [22]. The IR poses a major challenge to the efficacy of virotherapy
in eradicating viruses and CC. Clinical studies have revealed that the virus can activate the immune system against
the tumor, eradicating it [23,24]. Thus, we have formulated a mathematical model that explains the dynamics of

9
Fig. 2. Time variation in the relative population sizes with initial values. Here parameter values are: r = 0.43, a = 7/10×10−9 , b =
2 × 10−8 , k = 1.4, d = 0.1, e = 0.0009, c = 10−8 , m = 1.6 × 10−9 , and n = 2

Fig. 3. Time variation in the relative population sizes with initial values. Here parameter values are: r = 0.45, a = 7/10 ×
10−9 , b = 2 × 10−8 , k = 2, d = 0.23, e = 0.999, c = 10−4 , m = 0.0005, and n = 1

Fig. 4. Time variation in the relative population sizes with initial values. Here parameter values are: k = 3.9, d = 0.29, e =
0.0009, m = 2.1932, and n = 0.2

10
Fig. 5. Time variation in the relative population sizes with initial values. Here parameter values are: r = 0.4, a = 7/10 × 10−9 , b =
2 × 10−8 , k = 2, d = 0.1, e = 0.315, c = 10−8 , m = 0.0005, and n = 0.0001

Fig. 6. Time variation in the cancer cell of model with varying burst size (e) with initial values (T = 0.5, V = 0.01, and I = 0.01)

Fig. 7. Time variation in the cancer cell of model with varying values of m with initial values (T = 0.5, V = 0.01, and I = 0.01)

11
virotherapy, including the impact of IR as a predator. In this paper, we numerically analyzed a model that incorporates
the impact of IR on CC. We varied important parameters to identify optimal values for better outcomes in CC.
In section 3, we have done proof LAS for all the critical points. The model contained LAS for dead critical point when
the natural death rate of CC is higher than the logistic growth rate of CC or (r < d), for cancer invasion equilibrium
m(r − d)
point when < 1, for immune response free equilibrium point when mT2 < n and e > aV2 , for virus free
nkr
equilibrium point when we apply RHC, and for co-existing critical point when we apply RHC, and f > (r + mT4 ).
The model contained GAS for dead equilibrium point when r < d , and for cancer invasion equilibrium point when
fr
1− < 0 and R0 < 1.
a
In numerical simulation, section 4, Fig. 1 shows the stability of the dead critical state. Fig. 2 represents the stability
of cancer invasion critical point. Fig. 3 shows the stability of immune response free critical point. Fig. 4 shows the
stability virus free critical point, and Fig. 5 shows the stability co-existing critical point. From our analysis, we have
found that the success of virotherapy in controlling tumor size is directly tied to the viral infection rate. In other
words, for virotherapy to be effective, the viral infection rate needs to be high. High viral infection rates are crucial
for virotherapy to effectively control tumor size. Figure 6 shows the time variation in cancer cell for different values
of burst sizes. The figure shows that the largest burst size reduces tumor size. Furthermore, we have found that a high
stimulation of IR significantly diminishes the tumor size. Figure 7 shows the time variation in cancer cell for different
values of stimulation immune response. The figure shows that a high stimulation immune response rate reduces tumor
size.

6 Conclusion

We devolved a mathematical model of reduction of cancer cell using virotherapy and immune response. Then we
verified the concept by presenting positive and bounded solutions. Moreover, the model’s stability was analyzed using
nonlinear systems theory. We obtained five critical points, namely: dead critical point, cancer invasion critical point,
immune response-free critical point, virus-free critical point, and co-existing critical point. These points represent dif-
ferent outcomes of the treatment. E0∗ represents success in treatment, while E1∗ , indicates failure in treatment. E2∗ to
virus dominance, E3∗ to immune dominance, and E4∗ to the coexistence of all state variables. We have determined the
local and global stability of critical points. We have demonstrated the stability using numerical simulations. Addition-
ally, we performed a parameter analysis on two parameters: burst size (e) and the stimulation rate of IR by CC (m).
We found that, high immune response rate and a large burst size reduce tumor cells. We can use virotherapy with other
therapies to eradicate cancer cells in the model successfully for future work.
Author Contributions:
Ambika: Conceptualization; Investigation; Methodology; Validation; original draft. Arvind Kumar Sinha: Conceptu-
alization; Formal analysis; Investigation; Supervision; Validation; Visualization; Writing - review and amp; editing.
Funding: Not applicable. Data availability statement: Not applicable to this article as no data sets were generated or
analyzed during the current study. Disclosure statements:The authors declare that there is no conflict of interest.

References
1. Cancer. 2022 Feb 3; 2024 Jan 22;. https://www.who.int/news-room/fact-sheets/detail/cancer.

12
 2. Marelli G, Howells A, Lemoine NR, Wang Y. Oncolytic viral therapy and the immune system: a double-edged sword against
cancer. Frontiers in immunology. 2018 Apr 26;9:866.
3. Andtbacka RH, Kaufman HL, Collichio F, Amatruda T, Senzer N, Chesney J, et al. Talimogene laherparepvec improves durable
response rate in patients with advanced melanoma. Journal of clinical oncology. 2015 Sep 1;33(25):2780-8.
4. Immune Response. 2022 May 23; 2024 Jan 22;. https://microbenotes.com/immune-response/.
5. Wodarz D. Gene therapy for killing p53-negative cancer cells: use of replicating versus nonreplicating agents. Human gene
therapy. 2003 Jan 15;14(2):153-9.
6. Bajzer Ž, Carr T, Josić K, Russell SJ, Dingli D. Modeling of cancer virotherapy with recombinant measles viruses. Journal of
theoretical Biology. 2008 May 7;252(1):109-22.
7. Tian JP. The replicability of oncolytic virus: defining conditions in tumor virotherapy. Mathematical Biosciences & Engineer-
ing. 2011 May 31;8(3):841-60.
8. Wodarz D. Viruses as antitumor weapons: defining conditions for tumor remission. Cancer research. 2001 Apr 15;61(8):3501-
7.
9. Ashyani A, Mohammadinejad H, RabieiMotlagh O. Hopf bifurcation analysis in a system for cancer virotherapy with effect
of the immune system. Jordan J Math Stat. 2016 Jun 1;9:93-115.
10. Phan TA, Tian JP, et al. The role of the innate immune system in oncolytic virotherapy. Computational and mathematical
methods in medicine. 2017;2017(1):6587258.
11. Al-Tuwairqi SM, Al-Johani NO, Simbawa EA. Modeling dynamics of cancer virotherapy with immune response. Advances
in Difference Equations. 2020 Dec 2020:1-26.
12. Al-Johani N, Simbawa E, Al-Tuwairqi S. Modeling the spatiotemporal dynamics of virotherapy and immune response as a
treatment for cancer. Commun Math Biol Neurosci. 2019 Oct 27;2019:Article-ID.
13. Simbawa E, Al-Johani N, Al-Tuwairqi S, et al. Modeling the spatiotemporal dynamics of oncolytic viruses and radiotherapy
as a treatment for cancer. Computational and Mathematical Methods in Medicine. 2020;2020(1):3642654.
14. Fujie K, Ito A, Winkler M, Yokota T. Stabilization in a chemotaxis model for tumor invasion. Discrete Contin Dyn Syst. 2016
Jan 1;36(1):151-69.
15. Viglialoro G, Woolley TE. Boundedness in a parabolic-elliptic chemotaxis system with nonlinear diffusion and sensitivity and
logistic source. Mathematical Methods in the Applied Sciences. 2018 Mar 30;41(5):1809-24.
16. Li T, Pintus N, Viglialoro G. Properties of solutions to porous medium problems with different sources and boundary condi-
tions. Zeitschrift für angewandte Mathematik und Physik. 2019 Jun;70:1-18.
17. Edelstein-Keshet L. Mathematical models in biology. Society for Industrial and Applied Mathematics; 2005.
18. Martcheva M. An introduction to mathematical epidemiology. New York:Springer. vol. 61; 2015.
19. Friedman A, Tian JP, Fulci G, Chiocca EA, Wang J. Glioma virotherapy: effects of innate immune suppression and increased
viral replication capacity. Cancer research. 2006 Feb 15;66(4):2314-9.
20. Titze MI, Frank J, Ehrhardt M, Smola S, Graf N, Lehr T. A generic viral dynamic model to systematically characterize
the interaction between oncolytic virus kinetics and tumor growth. European Journal of Pharmaceutical Sciences. 2017 Jan
15;97:38-46.
21. Tao Y, Guo Q. A free boundary problem modelling cancer radiovirotherapy. Mathematical Models and Methods in Applied
Sciences. 2007 Aug;17(08):1241-59.
22. Virotherapy. 2024 Jan 22;. https://en.wikipedia.org/wiki/Virotherapy.
23. Ribas A, Dummer R, Puzanov I, VanderWalde A, Andtbacka RH, Michielin O, et al. Oncolytic virotherapy promotes intratu-
moral T cell infiltration and improves anti-PD-1 immunotherapy. Cell. 2017 Sep 7;170(6):1109-19.
24. Kim M, Nitschké M, Sennino B, Murer P, Schriver BJ, Bell A, et al. Amplification of oncolytic vaccinia virus widespread
tumor cell killing by sunitinib through multiple mechanisms. Cancer research. 2018 Feb 15;78(4):922-37.

13
 Declaration of original manuscript

Dated 26 July 2024

Anthony J. McGoron,
Editor-In-Chief,
Critical Reviews™ in Biomedical Engineering. I have the pleasure of sending you the manuscript entitled “Mathe-
matical model of cancer treatment with virotherapy and immune system,” authored by Ambika and Dr. Arvind Kumar
Sinha, to be considered for publication as an original research article in your prestigious and esteemed journal Critical
Reviews™ in Biomedical Engineering. The Paper contains original research and has not been submitted/published
earlier in any journal and is not being considered for publication elsewhere. All authors have seen and approved
the manuscript and have contributed significantly to the paper. The article has also been spell-checked and carefully
checked for grammar by Grammarly.
Abstract
Cancer affects people worldwide. Although traditional treatments like chemotherapy and radiation are effective, they
have side effects and harm healthy cells. Virotherapy is a treatment that specifically targets cancer cells, leaving
healthy cells unharmed. This paper presents a mathematical model to explore the dynamic interplay among cancer, vi-
rotherapy, and the immune response. The model considers virus replication within the tumor, immune cell infiltration,
and the tumor-immune microenvironment. We verify the model by showing that the solutions are non-negative and
bounded. Our model contains five equilibrium points: dead equilibrium point, cancer invasion equilibrium point, im-
mune response-free equilibrium point, virus-free equilibrium point, and coexisting equilibrium point. We find the local
stability and global stability of the equilibrium points. We simulate numerically to demonstrate our analytical results.
A high stimulation immune response rate and a large burst size reduce tumor size. The paper presents the importance
of high immune system and virotherapy on cancer patients. The model helps to oncologist and immunologists to find
the behavior of virotherapy and immune response to control the proliferation of different kinds of tumors.
Ethical Procedure
• The research meets all applicable standards with regard to the ethics of experimentation and research integrity, and
the following are being certified/declared true.
• As an author and along with co-authors of the concerned field, the paper has been submitted with full responsibility,
following due ethical procedure, and there is no duplicate publication, fraud, plagiarism, or concerns about animal or
human experimentation.
1∗
Corresponding Author:
Ambika,
Department of Mathematics
National Institute of Technology Raipur
Raipur – 492010, Chhattisgarh (India)
E-mail: ambika11sahu@gmail.com
Phone No.: +91-9098141269

14
 Highlights: Mathematical model of cancer treatment with virotherapy and immune system

• Exploration of nonlinear trends of virotherapy and immune response.

• Obtain stability on immune response-free, virus-free and present equilibriums.
• A continuous model using ordinary differential equations.
• Obtain parameter analysis.
• Reduces the tumor size with burst size and stimulation rate of immune response.

15
 Conflict of interest statement

• None of the authors of this paper has a financial or personal relationship with other people or organizations that could
inappropriately influence or bias the content of the paper.
• It is to specifically state that “No Competing interests are at stake, and there is No Conflict of Interest” with other
people or organizations that could inappropriately influence or bias the content of the paper.
1∗
Corresponding Author:
Ambika,
Department of Mathematics
National Institute of Technology Raipur
Raipur – 492010, Chhattisgarh (India)
E-mail: ambika11sahu@gmail.com
Phone No.: +91-9098141269

16