Optimal Immunotherapy of Oncolytic Viruses and Adopted Cell

Transfer in Cancer Treatment

G.V.R.K. VITHANAGE, SOPHIA R-J JANG

Department of Mathematics and Statistics

Texas Tech University

1108 Memorial Circle, Lubbock, TX 79409-1042

USA

Abstract: We investigate therapeutic effects of monotherapy of oncolytic viruses, of adopted cell transfer, as well

as the two combined therapies over a short time treatment period by applying optimal control techniques. The

goal is to minimize the number of susceptible tumor cells and the costs associated with the therapy over the

treatment period. We verify that there exists an optimal control pair and derive the necessary conditions. The

optimality system is solved numerically to provide optimal protocols under different scenarios with respect to

initial tumor sizes and parameter values. Although the two types of therapy do not work synergistically when

the viral killing rate by immune cells is large, a small anti-viral killing can improve therapy success of either

monotherapy of oncolytic viruses or combined therapy of oncolytic viruses and adopted T cell transfer. This

finding can be accomplished either by manipulating certain genes of viruses via genetic engineering or by chemical

modification of viral coat proteins to avoid detection by the immune cells.

Key-Words: Tumor, Immune system, Oncolytic virus therapy, adoptive cell transfer, optimal control theory,

Pontryagin’s maximum principle

Received: May 27, 2021. Revised: March 18, 2022. Accepted: April 19, 2022. Published: June 7, 2022.

1 Introduction

Cancer immunotherapy is a promising strategy to

treat cancer of various types. It consists of activat-

ing and harnessing the immune system to fight can-

cers [7, 31]. The therapy takes on several differ-

ent approaches such as immune checkpoint inhibitors,

adoptive cell transfer, monoclonal antibodies, cancer

vaccines, and oncolytic virus therapy (OVT) [7, 31].

Adoptive cell transfer (ACT) is implemented through

selecting and expanding patients’ own tumor-specific

T cells in vitro and then reinfusing back into pa-

tients to boost the patients’ own immune ability to

target cancer [16, 33]. Oncolytic viruses (OVs), on

the other hand, are genetically engineered or natu-

rally occurring viruses that selectively replicate in

and kill cancer cells without harming the normal

tissues [12, 31]. In recent years, an array of on-

colytic viruses have demonstrated anti-tumor effi-

cacy, including adenoviruses, herpes simplex viruses,

measles viruses, vesicular stomatitis virus and New-

castle disease virus [5, 14, 17, 31]. Compared with

conventional treatment strategies, immune system ac-

tivating agents produced by oncolytic viruses enable

the infected tumor cells to be localized and concen-

trated, reducing possible side effects. In 2015, the

USA Food and Drug Administration (FDA) approved

Talimogene Laherparepvec (T-VEC), an herpes sim-

plex virus, as the first oncolytic virus for the treatment

of advanced melanoma [1, 10]. A wide variety of OVs

are currently going through studies in phase I/II clin-

ical trials or in preclinical cancer models [1].

The interaction between tumor and tumor mi-

croenvironment is a very complicated process and this

activity is changing rapidly. Mathematical modeling

provides a valuable tool for understanding the com-

plex interaction among the many components of the

tumor microenvironment [8, 13]. Mathematical mod-

els can represent a natural phenomena by an equation

system. This process can allow a more efficient and

accurate analysis of the phenomenon. In addition to

clinical and experimental research, mathematical and

computational models can be very useful in studying

various mechanisms involved in cancer development

and can be used to illuminate the underlying dynamics

of therapy systems, which can lead to more optimal

treatment strategies.

The immune system interacts intimately with tu-

mors over the entire process of disease development

and progression. This complex cross talk between im-

munity and cancer cells can both inhibit and enhance

tumor growth [7, 15, 36]. Recently, Vithanage et al.

[37] proposed and investigated a mathematical model

of tumor-immune system interactions with oncolytic

viral therapy, wherein the immune cells can either be

simulated to proliferate or be suppressed to increase

their mortality. They concluded that reducing the vi-

ral killing rate by immune cells always increases the

effectiveness of the viral therapy. The reduction of

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

140

Volume 19, 2022

anti-viral killing rate can be obtained by either ma-

nipulating certain genes of viruses via genetic engi-

neering or by chemical modification of viral coat pro-

teins to evade recognition by the immune cells [37].

However, their conclusion is based on the long term

dynamical behavior of the tumor by analyzing the dy-

namical system of ordinary differential equations both

analytically and numerically.

The goal of this work is to study short term effects

of the therapy using optimal control theory. Several

researchers have applied optimal control techniques

to study OVT strategies [23, 28, 35]. For example,

Kim et al. [23] investigated a tumor-immune system

model with viral and immunotherapy in an optimal

control setting. The immune cells in their model do

not kill viruses and the infection rate is modeled by a

simple mass action. Malinzi et al. [28] applied opti-

mal control theory to study a model of tumor-immune

system interactions with viral therapy and chemother-

apy and they demonstrated numerically that viral ther-

apy can enhance chemotherapy. Su et al. [35], on

the other hand, studied a model with oncolytic viral

therapy along with MEK inhibitor and compared op-

timal control strategy on the dosage of MEK inhibitor

with the constant control strategy. Their simulations

showed that the optimal control has better control ef-

fect than constant control.

In this investigation, the immune cells can kill

viruses from the oncolytic therapy and the viral infec-

tion rate is modeled by a Michaelis-Menten kinetics.

Specifically, the model of no optimal control is based

on the work of Vithanage et al. [37] in which there

are three boundary equilibria whereas the number of

positive equilibria cannot be determined analytically.

We will devise a best immunotherapy using optimal

control theory techniques. In the following section,

the mathematical model in [37] is revisited. An op-

timal control problem is introduced in Section 3 and

the existence of an optimal control pair is verified. We

apply the Pontryagin’s maximum principle to derive

the necessary conditions and the optimality system in

the subsections. Section 4 presents numerical exam-

ples to illustrate monotherapy of OVT, of ACT, and a

combination of OVT and ACT. The paper ends with a

summary and conclusions in the final section.

2 A Mathematical Model

We first provide a brief description of the model stud-

ied in [37].

The tumor cells are classified into susceptible x

and infected y. The compartment of free viral parti-

cles is denoted by v, and zis the collection of immune

cells including both the innate and adaptive immune

cells. The units of cell populations are numbers of

cells and the unit of viral particles is given by pfu, the

plaque-forming unit, with day as the time unit.

The susceptible tumor population xgrows logis-

tically with intrinsic growth rate rand carrying ca-

pacity 1/bfor all tumor cells [26, 30, 34, 38, 40].

The mass action kinetics has been used by numer-

ous researchers [22, 32, 34, 38] to describe viral in-

fection of susceptible tumor cells. However, virus

spread is likely to be slower, limited by spatial con-

straints [39]. Since viruses released from one infected

cell cannot reach all susceptible tumor cells, the in-

fection rate must be a saturating function of suscep-

tible tumor cells [39]. As in [37], the infection pro-

cess is modeled by a Michaelis-Menten term with a

half saturation constant gand a maximum infection

rate β. How the immune cells detect cancer cells

and kill them is a complicated process [7]. It is as-

sumed that the tumor killing by immune cells follows

a mass action law [13, 25] and let kdenote the rate by

which susceptible tumor cells are killed by immune

cells. The corresponding killing rate of infected tu-

mor cells by immune cells is given by c. The infected

tumor cells have an extra death rate adue to infection

[30, 32, 34, 38].

Upon lyses, new progeny virions are released from

infected tumor cells and the viral burst size per in-

fected tumor cell is denoted by q[30, 32, 34, 38]. Im-

mune cells recognize viruses as foreign pathogens and

mount strong anti-viral responses which can eventu-

ally kill viruses from the OVT. We model this killing

by the law of mass action with rate γ[3]. The natural

death rates of viruses and immune cells are constant

and are denoted by δand d, respectively [34, 38]. As

in [13, 25], it is assumed that there is a constant sup-

ply of immune cells from the lymph nodes into the tu-

mor microenvironment at a rate s. When OVs infect

and destroy cancer cells, specific antigens are released

into the tumor microenvironment allowing immune

cells to be activated [14, 29]. We use a Michaelis-

Menten mechanism to model this activation from in-

fected tumor cells [22, 30, 34] with pdenoting the

maximum rate and hthe half saturation constant.

There are numerous mechanisms used by tumors

to counter effect immune cells [7, 19]. Therefore, im-

mune cells can either be stimulated to proliferate or be

impaired to reduce their growth by susceptible can-

cer cells. A combination of stimulation and suppres-

sion on growth between cancer and immune cells may

act simultaneously and may be modeled by the sim-

ple mass action at a rate m[13, 37]. The net growth

increases or decreases due to presence of cancer cells

can either be positive (m > 0), negative (m < 0) or

null (m= 0). Negative mindicating tumor cells ex-

ert significant effect on the immune system leading to

immunosuppression.

Based on the above descriptions, the model takes

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

141

Volume 19, 2022

the following form

x0=rx1−b(x+y)−βxv

g+x−kxz

y0=βxv

g+x−ay −cyz

v0=qay −δv −γvz

z0=s−dz +mxz +pyz

h+y

(1)

with the initial condition x(0) >0, y(0) ≥0, x(0)+

y(0) <1/b, v(0) ≥0, z(0) ≥0. All of the parame-

ters r, b, β, k, a, c, q, δ, γ, s, d, p, h are positive

except mwhich can be any real number. The negative

mindicating the net effect of tumor on immune cells

decreases their growth.

The authors in [37] provided several sufficient

conditions in terms of model parameters for which

the tumor can be eradicated for all sizes, that is,

lim

t→∞ x(t)=0for all 0< x(0) <1/b. In particular,

if m < 0and k > r(bd −m)

sb , then the tumor-free

equilibrium E0= (0,0,0, s/d)is globally asymptot-

ically stable [37, Theorem 3.3]. However, it may take

a very long time, longer than a patient’s life time, to

eliminate the tumor.

3 Optimal Control Problem

As elaborated above, the eradication results derived in

[37] are with respect to the long term dynamics of the

tumor, which may not be applicable in real life sce-

nario. The aim of this section is to apply optimal con-

trol techniques to understand short term effects of im-

munotherapy relative to tumor regression. The goal is

to minimize the number of susceptible tumor cells and

the costs associated with immunotherapies over the fi-

nite treatment period [0, T ], where T > 0is fixed.

Let s1≥0and s2≥0be the strengths of oncolytic

viral therapy (OVT) and immunotherapy of adoptive

cell transfer (ACT) respectively with s1+s2>0.

Denote u1(t)and u2(t)the controls for the OVT and

ACT, respectively. In particular, the units of s1and s2

are pfu day−1and cell day−1respectively, whereas u1

and u2are dimensionless. The state equations take the

following form

x0=rx1−b(x+y)−βxv

g+x−kxz

y0=βxv

g+x−ay −cyz

v0=qay −δv −γvz +s1u1(t)

z0=s−dz +mxz +pyz

h+y+s2u2(t)

(2)

and the initial condition are given by x(0) >

0, y(0) ≥0, x(0)+y(0) <1/b, v(0) ≥0, z(0) ≥0.

Since the goal is to minimize susceptible tumor size

as well as the costs of implementing immunotherapies

over the treatment period [0, T ], the objective func-

tional is written as

J(u1, u2) = ZT

(x(t) + c1

2(u1(t))2+c2

2(u2(t))2)dt,

(3)

where (u1, u2)belonging to the class

U={(u1, u2) : ui(t)is piecewise continuous with

0≤ui(t)≤1on [0, T ], i = 1,2}.(4)

The parameters c1≥0and c2≥0are the weighted

constants used to balance the contributions between

the two types of treatment. We assume ci>0if si>

0. The optimal control problem consists of

min

(u1,u2)∈UJ(u1, u2)(5)

subject to the state equations (2).

Once the optimal control problem is specified, we

proceed to discuss the existence of optimal controls,

their characterizations and the optimality system in

the following subsections.

3.1 Existence of optimal control pair

Classical theory of optimal control [11] can be applied

directly to study the problem formulated in (2)–(5).

We first verify the existence of an optimal control pair.

Theorem 3.1. There exists an optimal control pair for

the problem (2)–(5).

Proof. It is enough to show that the following condi-

tions given in Corollary 4.1 of [11] are satisfied.

(a) The set of all initial conditions with a control pair

(u1, u2)∈Ufor which the state equations being

satisfied is nonempty.

(b) Uis closed and convex.

is continuous, bounded above by the sum of the

control and the state, and can be written as a lin-

ear function of u1(t),u2(t)with coefficients de-

pending on time and the state.

(d) The integrand of J(u1, u2)is convex in Uand

is bounded below by −k2+k1|(u1, u2)|ηwith

k1>0and η > 1.

Clearly for each fixed initial condition and control

pair, (1) has a unique solution on [0, T ]and (a) is

satisfied. Moreover, as x0|x=0 = 0,y0|y=0 ≥0,

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

142

Volume 19, 2022

v0|v=0 ≥0and z0|z=0 >0, solutions remain non-

negative on [0, T ]. It is obvious that (b) is true and

(d) is satisfied with η= 2. To verify (c), notice

x(t)>0, y(t)≥0, v(t)≥0, z(t)≥0and

0< x(t) + y(t)<1/bfor all t∈[0, T ]. Thus

x0≤rx,y0≤βv −ay,v0≤qay −δv +s1u1(t),

and z0≤s−dz +pz +s2u2(t)for m < 0while

z0≤s−dz +mz/b+pz +s2u2(t)if m≥0, i.e.,













≤M











+





s1u1(t)

s+s2u2(t)







where

M=





r0 0 0

0−a β 0

0qa −δ0

0 0 0 −d+m/b+p







if m≥0

and

M=





r0 0 0

0−a β 0

0qa −δ0

0 0 0 −d+p







if m < 0.

Let X= (x, y, v, z)tr, the transpose of (x, y, v, z).

Then

||dX

dt || ≤ ||M|| · || 











|| +|| 





s1u1

s+s2u2







||.

Thus (c) is verified and there exists an optimal control

pair for the control problem (2)–(5) by [11, pages 68-

69].

3.2 The adjoint system and characterization

of control pair

We apply the Pontryagin’s maximum principle to de-

rive necessary conditions [27]. Let (λ1, λ2, λ3, λ4)

denote the adjoint vector. The Hamiltonian of the op-

timal control problem (2)–(5) is

H(x, y, v, z, u1, u2, λ1, λ2, λ3, λ4)

=x+c1

2u2

1+c2

2u2

+λ1rx(1 −b(x+y)) −βxv

x+g−kxz

+λ2βxv

x+g−ay −cyz

+λ3qay −δv −γvz +s1u1

+λ4s−dz +mxz +pyz

h+y+s2u2

(6)

where the adjoint variables satisfy

λ0

1=−∂H

∂x , λ0

2=−∂H

∂y ,

λ0

3=−∂H

∂v , λ0

4=−∂H

∂z ,

with the transversality conditions

λi(T) = 0 for 1≤i≤4.

Setting ∂H

∂ui

= 0,i= 1,2, we obtain u1=

−λ3s1

and u2=−λ4s2

. Since the controls u1and

u2are bounded, 0≤u1, u2≤1, the characterization

of the optimal control pair is therefore given by

u∗

1(t) = min 1,max{0,−s1λ3

}

u∗

2(t) = min 1,max{0,−s2λ4

}(7)

provided ci>0,i= 1,2. The above discussion is

summarized as follows.

Proposition 3.1. Given an optimal control pair

(u∗

1, u∗

2)and solutions of the corresponding state

equations (2), there exist adjoint variables λi,1≤

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

143

Volume 19, 2022

i≤4, satisfying

λ0

1=−1−βgv

(g+x)2λ2−mzλ4

−r(1 −2bx −by)βgv

(g+x)2−kzλ1

λ0

2=λ1rbs + (a+cz)λ2−qaλ3−phz

(h+y)2λ4

λ0

3=βx

g+xλ1−βx

g+xλ2+ (δ+γz)λ3(8)

λ0

4=kxλ1+cyλ2+γvλ3

−(−d+mx +py

h+y)λ4

λi(T) = 0,1≤i≤4.

Moreover, u∗

1and u∗

2are represented by (7).

The optimality system consisting of the state and

adjoint equations is given as

x0=rx1−b(x+y)−βxv

g+x−kxz

y0=βxv

g+x−ay −cyz

v0=qay −δv −γvz

+s1min 1,max{0,−s1λ3

}

z0=s−dz +mxz +pyz

h+y

+s2min 1,max{0,−s2λ4

}

λ0

1=−1−λ1r(1 −2bx −by)−βgv

(g+x)2−kz

−λ2

βgv

(g+x)2−mzλ4(9)

λ0

2=rbxλ1+ (a+cz)λ2−qaλ3−phz

(h+y)2λ4

λ0

3=λ1

βx

g+x−λ2

βx

g+x+λ3(δ+γz)

λ0

4=hxλ1+cyλ2+γvλ3

−λ4(−d+mx +py

h+y)

with boundary conditions x(0) >0, y(0) ≥

0, v(0) ≥0, z(0) >0, x(0) + y(0) <1/b, λi(T) =

0,1≤i≤4. The optimality system (9) yields a

two-point boundary value problem. It can be verified

as in [2, 20] that the solution of (9) is unique if T > 0

is small.

4 Numerical Investigations

We apply the numerical technique of backward-

forward sweep method described in [27] combined

with the fourth-order Runge-Kutta scheme to numer-

ically solve the optimality system (9).

Before numerical explorations, a plausible range

of parameter values and their sources are provided in

Table 1 with the baseline values given in equation

(10). There is a wide range of tumor growth rates

in the literature. For example, r= 0.18 in [25],

r= 0.514 in [6, 13], and r∈(0.69,0.97) are adopted

in [9]. In the review paper [9], tumor growth rates of

0.23,0.43 and 1.636 are simulated in the numerical

examples to demonstrate model outcomes. Motivated

by the work in [21], the tumor growth rate is estimated

using the doubling time of about 48-60 hours of the

BxPC-3 cell line [4] by assuming an initial exponen-

tial growth phase. These parameter values were also

adopted in the numerical investigations in [37].

Parameter Value Reference

r0.2773–0.3466 day−1[37]

b1.02×10−9cell−1[6]

β6×10−12–0.862 cell pfu−1day−1[9]

k10−5–10−3cell−1day−1[13]

a1.333–2.6667 day−1[13]

c0.0096–4.8 cell−1day−1[13]

q10–1350 pfu day−1[13]

δ0.024–24 day−1[34]

γ0.024–48 cell−1day−1[34]

m-1–1.5×10−9cell−1day−1[13]

p2.4×10−4–2.5008 day−1[34]

h20–5×104cell [34]

s5×103cell day−1[13]

d2 day−1[13]

g40–105cell [30]

s1varied pfu day−1

s2varied cell day−1

u10–1 dimensionless [27]

u20–1 dimensionless [27]

Table 1: Parameters values and their sources

The baseline parameter values are given by

r= 0.346, b = 1.02 ×10−9, a = 1.333, c = 1.8,

q= 100, d = 2, δ = 1.83, p = 2.4×10−4,(10)

h= 5 ×104, g = 105, s = 5 ×103,

and we vary the parameter values of β,γ,mand k.

Unless otherwise stated, c1= 10 and c2= 10 are

used when the corresponding therapy is applied. The

values of c1and c2are chosen arbitrary. However,

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

144

Volume 19, 2022

if we vary these values we obtain similar simulation

results.

In [24], 5×106human BxPC-3 cells of pancre-

atic ductal adenocarcinoma were injected into mice

subcutaneously. When the tumor had grown to a di-

ameter of 5–7 mm, two different oncolytic viruses at

a MOL of 108pfus were injected subcutaneously to

different groups of mice. Therefore, the dosage s1of

OVT is varied around 108pfus. In [18], 105to 2×106

number of T cells were infused to the experimental

mouse. In our numerical examples, the values of s1

and s2are hypothetical but are close to the values of

above experiments.

We first consider the monotherapy of OVT, of

ACT, and then the combined immunotherapy of OVT

and ACT in the following subsections.

4.1 Monotherapy of oncolytic viruses

We explore the situation when only the OVT is ap-

plied by letting c2=s2= 0.

First, consider the strength s1of OVT as s1= 108

with a fixed treatment period of 100 days, i.e., T=

100. Unless otherwise stated, the initial condition is

chosen as

(x(0), y(0), v(0), z(0)) = (2 ×107,0,0,300),

which is the larger tumor size given in the numerical

examples of [37]. The simulation results are provided

in Figures 1–4, where the susceptible tumor sizes of

no therapy denoted by x(t)and of optimal OVT repre-

sented by x∗(t)are plotted in the top row. The bottom

row plots the corresponding optimal controls u∗

i(t).

Figure 1 shows that the OVT has to be performed dur-

ing the whole treatment period in order to control the

tumor. Notice that the susceptible tumor is increas-

ing in Figure 1(b) and 1(c) when there is no therapy,

and it is either stabilized in Figure 1(b) or slowly de-

creases in Figure 1(c) during the whole treatment pe-

riod. With the parameter values in Figure 1(d), the

tumor size is decreasing slowly when no OVT is ap-

plied. However, the OVT has to be applied for the

whole treatment period in order to reduce the suscep-

tible tumor size more effectively.

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

0 20 40 60 80 100

2.5 107

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

(a) (b)

0 20 40 60 80 100

2.5 107

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

0 20 40 60 80 100

2.2

2.4 107

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

Figure 1: s1= 108and β= 0.6with γ= 0.9in (a),

(b) and (d), and γ= 0.7in (c). (a) m=−10−6, k =

10−5. (b)–(c) m= 10−9, k = 1.33 ×10−4. (d)

m= 1.5×10−9,k= 1.33 ×10−4.

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

0 20 40 60

107

x*(t)

x(t)

0 20 40 60

days

0.5

(a) (b)

0 10 20 30 40

107

x*(t)

x(t)

0 10 20 30 40

days

0.5

0 5 10 15 20

10 108

x*(t)

x(t)

0 5 10 15 20

days

0.5

Figure 2: Parameters s1= 108and γ= 0.05 are

fixed. (a) m=−10−6,β= 0.6,k= 1.33 ×10−4.

(b) m= 10−9,β= 0.6,k= 1.33 ×10−4with

T= 65. (c) m= 10−9,β= 0.8,k= 1.33 ×10−4.

(d) m=−10−6,β= 0.8,k= 10−3.

Next, s1= 108with γ= 0.05 are simulated and

the results are provided in Figure 2. The viral killing

rate by immune cells is smaller than in Figure 1 and

we can see that the OVT is more effective. Particu-

larly, the treatment period can be shortened to 65 days

in Figure 2(c) as the susceptible tumor cells are erad-

icated by day 60 using the OVT alone. Comparing

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

145

Volume 19, 2022

Figure 2(c) and (d), there is a change in tumor killing

rate kbut a huge difference between the two treatment

outcomes even though the tumor microenvironment is

more immunosuppressive in Figure 2(d).

In Figure 3 we choose the viral infection rate βas

0.8and let γ= 0.05 with various values of mand

k. When the tumor killing rate kis increased from

10−4to 1.33 ×10−4, the susceptible tumor is eradi-

cated by day 45 shown in Figure 3(a) and (b). If tumor

is more immunosuppressive but with a larger tumor

killing rate k= 10−3, the tumor can be eliminated by

day 10 as illustrated in Figure 3(d). On the contrary,

the tumor cannot be eradicated if the tumor killing rate

is smaller with k= 10−4in Figure 3(c).

0 20 40 60 80 100

2108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

0 10 20 30 40

107

x*(t)

x(t)

0 10 20 30 40

days

0.5

(a) (b)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

0 5 10 15 20

10 108

x*(t)

x(t)

0 5 10 15 20

days

0.5

Figure 3: s1= 108,β= 0.8and γ= 0.05. (a) m=

10−9,k= 10−4. (b) m= 10−9,k= 1.33 ×10−4

with T= 45. Notice that x∗(t)becomes extinct by

t= 45. (c) m=−10−6,k= 10−4. (d) m=−10−6,

k= 10−3,T= 20.

In the next figure, Figure 4, we let β= 0.6,

m=−10−6but with a smaller initial tumor size

(107,0,0,300). In this scenario, the susceptible tu-

mor is oscillating over the treatment period with sizes

smaller than the tumor with no treatment. In addition,

the therapy needs not to be implemented for the whole

treatment period. There are short periods of times

when no therapy is required. We also simulated the

cases that correspond to Figure 4 with c1= 100. Sim-

ilar results are obtained and they are not presented.

(a) (b)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

Figure 4: s1= 106,β= 0.6,m=−10−6, and

(x(0), y(0), v(0), z(0)) = (107,0,0,300). (a) k=

10−5,γ= 0.15 . (b) k= 10−4,γ= 0.15. (c)

k= 10−5,γ= 0.25. (d) k= 10−4,γ= 0.25. The

susceptible tumor without treatment and with OVT

are given by x(t)and x∗(t)respectively.

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

(a) (b)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

Figure 5: s1= 106,m=−10−6,k= 10−5, and

(x(0), y(0), v(0), z(0)) = (107,0,0,300) are fixed.

(a) β= 0.58,γ= 0.15. (b) β= 0.61,γ= 0.15. (c)

β= 0.6,γ= 0.16. (d) β= 0.6,γ= 0.17. The sus-

ceptible tumor without treatment and with OVT are

given by x(t)and x∗(t)respectively

Figure 5 presents simulation results with a smaller

tumor size x(0) = 107under different scenarios. The

only difference between (a) and (b) is the viral infec-

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

146

Volume 19, 2022

tion rate, where βis 0.58 in (a) and 0.61 in (b), and

their effect is small with respect to the timing of ther-

apy as well as susceptible tumor size. In (c) and (d),

the viral infection rate βis the same but with differ-

ent viral killing rate γ. It is 0.16 in (c) and 0.17 in (d).

With a slight larger viral killing rate, the viral therapy

will need to apply for a longer period of time in order

to control the tumor.

4.2 Monotherapy of adopted cell transfer

In this subsection, monotherapy of adopted immune

cell transfer (ACT) is considered, i.e., c1=s1= 0.

If there are no virus and infected tumor cells initially,

v(0) = 0 = y(0), then since only ACT is applied, it is

sufficient to consider the xz-subsystem with ACT. In

particular, parameter values of βand γwill not affect

treatment outcome.

012345

10 107

x*(t)

x(t)

012345

days

0.5

012345

4108

x*(t)

x(t)

012345

days

0.5

(a) (b)

0 5 10 15

10 108

x*(t)

x(t)

0 5 10 15

days

0.5

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

days

0.5

Figure 6: (a) m=−10−6,k= 10−5,

s2= 108, and (x(0), y(0), v(0), z(0)) = (2 ×

107,0,0,300). (b) m=−10−6,k= 10−5,

s2= 108,(x(0), y(0), v(0), z(0)) = (108,0,0,300).

(x(0), y(0), v(0), z(0)) = (2 ×107,0,0,300).

(d) m= 10−9,k= 10−5,s2= 105,

(x(0), y(0), v(0), z(0)) = (108,0,0,300). The sus-

ceptible tumor without treatment and with OVT are

given by x(t)and x∗(t)respectively

The parameter values m=−10−6and k= 10−5

are fixed in Figure 6(a)-(b) and we vary s2and the ini-

tial susceptible tumor size. In this scenario, the crit-

ical immune cell size needed to eradicate susceptible

tumor cells derived in [37] is 3.3991×107. From Fig-

ure 6(a), we see that with an ACT strength of s2= 108

the x∗(t)can be eliminated within the first 5 days.

We increase the initial tumor size x(0) to 108in Fig-

ure 6(b) and obtain a similar conclusion. In Figure

6(c)-(d), the tumor microenvironment is less suppres-

sive with m= 10−9. The result shows that a smaller

tumor size can be eradicated by day 100 using only

s2= 105.

4.3 Combined OVT and ACT

We consider combined therapy of OVT and ACT in

this subsection. In Figure 7(a)-(b), β= 0.6,k=

10−5,m=−10−6,s1= 105and s2= 103with

initial condition (2 ×107,0,0,300) are fixed. The

anti-viral killing rate γis 0.9in Figure 7(a) and 0.09

in Figure 7(b). We see that if γis reduced, the ACT is

more valuable and effective. The OVT needs not to be

implemented for the whole treatment period and the

tumor can be reduced further with oscillations under

more frequent ACT. In Figure 7(c), β= 0.9,k=

10−5,m=−10−6,γ= 0.09,s1= 104and s2= 102

with initial condition (2 ×107,0,0,300). The virus

is more transmissible but with a smaller dose of OVT.

The OVT is applied for the whole treatment period in

order to control the tumor.

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

(a) (b)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

(c)

Figure 7: k= 10−5,m=−10−6with initial condi-

tion (2 ×107,0,0,300). Then s2cri = 3.3991 ×107.

(a) β= 0.6, γ = 0.9,s1= 105,s2= 103. (b)

β= 0.6,γ= 0.09,s1= 105,s2= 103. (c) β= 0.9,

γ= 0.09,s1= 104,s2= 102.

In Figure 8, β= 0.9,γ= 0.09,m=−10−6,s1=

104,s2= 102and initial condition (2×107,0,0,300)

are fixed. We vary the anti-tumor killing rate kwith

k= 8 ×10−4,7×10−4and 6×10−4in Figure 8(a),

(b) and (c), respectively. Since the anti-vrial killing

rate γ= 0.09 is small, ACT does not affect much of

the OVT efficacy. The OVT has to be applied for the

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

147

Volume 19, 2022

whole treatment period when kis small. The therapy

does not need to be employed for the whole treatment

period when kis larger.

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

(a) (b)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

(c)

Figure 8: β= 0.9,γ= 0.09,m=−10−6,s1= 104,

s2= 102with initial condition (2 ×107,0,0,300)

are fixed. (a) k= 8 ×10−4. (b) k= 7 ×10−4.(c)

k= 6 ×10−4.

Comparing Figures 7 and 8, notice that when γ=

0.9is large in Figure 7(a), the ACT is detrimental to

the effectiveness of OVT and thus ACT is only ap-

plied initially and the tumor is not reduced signifi-

cantly. When γgets smaller with γ= 0.09 in Fig-

ure 7(b), the ACT and OVT work somewhat syner-

gistically and the tumor can be controlled with oscil-

lations. The tumor can be reduced further if β= 0.9

given in plot (c) of Figure 7.

Finally, k= 10−5,m=−10−6,s1= 105,s2=

103with initial condition (2×107,0,0,300) are fixed

while varying βand γin Figure 9. Comparing plots

(a) and (b), as γis increased, the ACT is detrimental

to the effectiveness of OVT and thus the ACT is only

applied in the beginning of the treatment period. As

γis increased further but with a larger infection rate

β, the ACT can be applied intermittently as shown in

plot (c).

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

(a) (b)

0 20 40 60 80 100

10 108

x*(t)

x(t)

0 20 40 60 80 100

0.5

*(t)

0 20 40 60 80 100

days

0.5

*(t)

(c)

Figure 9: k= 10−5,m=−10−6,s1= 105,s2=

103with initial condition (2×107,0,0,300) are fixed

for all the simulations. (a) β= 0.5,γ= 0.06. (b)

β= 0.5,γ= 0.07. (c) β= 0.8,γ= 0.1.

5 Summary and Conclusion

In this study, we apply optimal control techniques to

investigate effectiveness of the oncolytic viral therapy

and adopted T cell transfer over a short time treatment

period. Specifically, the model of no control is based

on an earlier study given in [37] where a mathematical

model of the tumor-immune system interactions with

OVT was studied. The authors in [37] derived several

sufficient conditions in terms of model parameters for

which the susceptible tumor cells can be eradicated

for all sizes. However, these results are based on the

asymptotic dynamics of the tumor which may not be

applicable in real life scenarios.

The goal of the control is to minimize the num-

ber of susceptible tumor cells along with the costs or

the tolerances associated with the therapy. We do not

minimize the number of infected tumor cells since

these cells do not proliferate. As immune cells can

kill viruses from the OVT, the two types of therapy do

not provide synergistic effects. This is especially true

when the viral killing rate γis large. When the dosage

s2of ACT is larger than the critical size derived in

[37], the susceptible tumor cells can be eradicated by

the first 20 days as shown in Figure 6. In general, the

monotherapy of OVT is more effective if the viral in-

fection rate βis larger while the viral killing rate γ

is smaller even when the tumor microenviroment is

very immunosuppressive. See Figures 1–3. If the ini-

tial tumor size is smaller, that is, x(0) = 107instead

of x(0) = 2 ×107, the OVT needs not to be imple-

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

148

Volume 19, 2022

mented for the whole treatment period and the suscep-

tible tumor exhibits dormant and elapse phenomenon

shown in Figures 4 and 5. For the combined therapy

of OVT and ACT, the ACT is not useful on control-

ling the tumor if anti-viral killing rate γis large as

shown in Figures 7–9. From these simulations, we

conclude that sole OVT is more effective than the sole

ACT if s2is below the critical value. However, if s2is

above the threshold, then the susceptible tumor cells

can be eradicated less than 100 days while the results

obtained in [37] are with respect to the long term dy-

namics of the susceptible tumor cells.

References:

[1] Apolonio J, de Souza Gonçalves V, Santos ML,

Luz M, Souza JV, Pinheiro SL, de Souza WR,

Loureiro M, de Melo FF. Oncolytic virus ther-

apy in cancer: A current review, World J Virol.,

210(5), 2021, 229-255.

[2] Burden, T., Ernstberger, J., Fister, K., Optimal

control applied to immunotherapy, Dis. Cont.

Dyn. Sys. Ser. B, 4, 2004,135-146.

[3] Choudhury B, Nasipuri B, Efficient virotherapy

of cancer in the presence of immune response, Int.

J. Dynam. Control,2, 2014, 314-325.

[4] Deer EL, et al., Phenotype and genotype of pan-

creatic cancer cell lines, Pancreas,39(4),2010,

425-435.

[5] de Matos AL, Lina S. Franco LS, McFadden G,

Oncolytic viruses and the immune system: The

dynamic duo, Mol. Ther. Methods Clin. Dev.,

17,2020, 349-358.

[6] de Pillis L, Radunskaya A, Wiseman C, A vali-

dated mathematical model of cell-mediated im-

mune response to tumor growth. Cancer Res.,

65(17),2005,7950-7958.

[7] Dong H, Markovic SN, The Basics of Cancer Im-

munotherapy, Springer 2018.

[8] Eftimie R, et al., Interaction between the immune

system and cancer: a brief review of non-spatial

mathematical models, Bull. Math. Biol.,73, 2011,

2-32.

[9] Eftimie R, Eftimie G, Tumour-associated

macrophages and oncolytic virotherapies: a

mathematical investigation into a complex

dynamics, Lett. Biomath. 5, 2018, 6-35.

[10] Ferrucci P, Pala L, Conforti F, Emilia Cocoroc-

chio E, Talimogene Laherparepvec (T-VEC):

An intralesional cancer immunotherapy for ad-

vanced melanoma, Cancers,13, 2021, 1383.

https://doi.org/10.3390/cancers13061383.

[11] Fleming, W., Rishel, R., Deterministic and

Stochastic Optimal Control, Springer, New York,

1975.

[12] Fukuhara H, Ino Y, Todo T, Oncolytic virus ther-

apy: A new era of cancer treatment at dawn. Can-

cer Sci.,107, 2016, 1373-1379.

[13] Garcia V, Bonhoeffer S, Fu F, Cancer-induced

immunosuppression can enable effectiveness of

immunotherapy through bistability generation: A

mathematical and computational examination, J.

Theor. Biol.,492, 2020, 110185.

[14] Gujar S, Pol JG, Kim Y, et al., Antitumor ben-

efits of antiviral immunity: An underappreci-

ated aspect of oncolytic virotherapies, Trends Im-

munol.,39, 2018, 209-221.

[15] Gun S, et al., Targeting immune cells for can-

cer therapy, Redox Biol., doi: 10.1016/j.re-

dox.2019.101174.

[16] Haddad D, Genetically engineered vaccinia

viruses as agents for cancer, treatment, imaging,

and transgene delivery, Front. Oncol.,7, 2017, 1-

12.

[17] Hale DF, Vreeland TJ, Peoples GE, Arming the

immune system through vaccination to prevent

cancer recurrence, Am. Soc. Clin. Oncol. Educ.

Book.,35, 2016, e159-e167.

[18] Hanada K, et al., An effective mouse model

for adoptive cancer immunotherapy targeting

neoantigens, JCI Insight, JCI Insight. 2019(10),

e124405.

[19] Hanahan D, Weinberg RA, Hallmarks of cancer:

the next generation, Cell,144(5), 2011, 646-674.

[20] Hu X, Jang R-J S, Optimal treatments in cancer

immunotherapy involving CD4+T cells, WSEAS

Trans. Biol. Biomed.,15, 2018, 48-67.

[21] Hu X, Ke G, Jang R-J S, Modeling pancre-

atic cancer dynamics with immunotherapy, Bull.

Math. Biol.,81, 2019, 1885-1915.

[22] Jang R-J S, Wei H-C, On a mathematical model

of tumor-immune system interactions with an on-

colytic virus therapy, Discrete Contin. Dyn. Syst.

Ser. B,27(6), 2022, 3261-3295.

[23] Kim KS, Kim S, Jung I, Hopf bifurcation analy-

sis and optimal control of Treatment in a delayed

oncolytic virus dynamics, Math. Comput. Simul.,

149, 2018, 1-16.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

149

Volume 19, 2022

[24] Koujima T, Tazawa H, Ieda T, et al., Oncolytic

virus-mediated targeting of the ERK signaling

pathway inhibits invasive propensity in human

pancreatic cancer, Mol. Ther. Oncolytics,17,

2020, 107-117.

[25] Kuznetsov VA, Makalkin IA, Taylor MA, Perel-

son AS, Nonlinear dynamics of immunogenic tu-

mors: parameter estimation and global bifurca-

tion analysis, Bull. Math. Biol.,56(2), 1994, 295-

321.

[26] Laird AK, Dynamics of tumor growth, Br. J.

Cancer,18, 1964, 490-502.

[27] Lenhart, L., Workman, JT., Optimal Control Ap-

plied to Biological Models, Chapman & Hall:

New York, 2007.

[28] Malinzi J, et al., Enhancement of chemotherapy

using oncolytic virotherapy: Mathematical and

optimal control analysis, Math. Biosci. Eng.,15,

2018, 1435-1463.

[29] Magen A, Nie J, Ciucci T, et al., Single-cell pro-

filing defines transcriptomic signatures specific

to tumor-reactive versus virus-responsive CD4+

T cells, Cell Reports,29, 2019, 3019-3032.

[30] Mahasa KJ, Eladdadi A, de Pillis L, Ouifki R,

Oncolytic potency and reduced virus tumorspeci-

ficity in oncolytic virotherapy. A mathematical

modelling approach, PLoS ONE,12(9), 2017, 1-

25.

[31] Marelli G, Howells A, Lemoine NR, Wang Y,

Oncolytic viral therapy and the immune system:

A double-edged sword against cancer, Front. Im-

munol.,9, 2018, 1-9.

[32] Okamoto K, Amarasekare P, Petty I, Model-

ing oncolytic virotherapy: Is complete tumor–

tropism too much of a good thing? J. Theor. Biol.,

358, 2014, 166-178.

[33] Rosenberg SA, et al., Adoptive cell transfer: a

clinical path to effective cancer immunotherapy,

Nat Rev Cancer,4(8), 2000, 10.1038/nrc2355.

[34] Storey KM, Lawler SE, Jackson TL, Model-

ing oncolytic viral therapy, immune checkpoint

inhibition, and the complex dynamics of innate

and adaptive immunity in glioblastoma treatment,

Front. Physiol.,11, 2020, 1-18.

[35] Su Y, Jia C, Chen Y, Optimal Control Model of

Tumor Treatment with Oncolytic Virus and MEK

Inhibitor, Biomed Res. Int., Volume 2016, 2016,

5621313.

[36] Vinay D et al., Immune evasion in cancer:

Mechanistic basis and therapeutic strategies,

Semin. Cancer Biol.,35, 2015, S185-S198.

[37] Vithanage R., Wei H-C, Jang S R-J., Bistabil-

ity in a model of tumor-immune system inter-

actions with an oncolytic viral therapy, Math.

Biosci. Eng.,19(2), 2022, 1559-1587.

[38] Wodarz D, Viruses as antitumor weapons, Can-

cer Res. 61, 2001, 3501-3507.

[39] Wodarz D, Komarova N, Towards predictive

computational models of oncolytic virus therapy:

Basis for experimental validation and model se-

lection, PLoS ONE,4, 2009, e4271.

[40] Wu JT, Byrne HM, Kirn DH, Wein LM, Mod-

eling and analysis of a virus that replicates selec-

tively in tumor cells, Bull. Math. Biol.,63, 2001,

731-768.

Contribution of individual authors to

the creation of a scientific article

(ghostwriting policy)

Conceptualization: Sophia Jang, Rohana Vithanage;

Formal analysis: Rohana Vithanage, Sophia Jang;

Numerical Simulations: Rohana Vithanage, Sophia

Jang;

Writing: Sophia Jang, Rohana Vithanage

Follow: www.wseas.org/multimedia/contributor-

role-instruction.pdf

Creative Commons Attribution

License 4.0 (Attribution 4.0

International , CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/li-

censes/by/4.0/deed.en_US

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.15

G. V. R. K. Vithanage, Sophia R-J Jang

E-ISSN: 2224-2902

150

Volume 19, 2022