Spiral waves and turbulence in mathematical models of oncolytic virotherapy induced by Hopf bifurcation

Oncolytic virotherapy is a promising targeted cancer treatment that employs viruses, which selectively infect tumor cells. Although its clinical efficiency has remained limited and it is often used in conjunction with other therapies, advances in genetic engineering have produced stronger and more selective viral strains, prompting continued interest in their dynamics. In particular, previous studies have noted that viruses with sufficiently high replication rates can induce oscillations reminiscent of predator-prey systems. Here, we extend this analysis to the spatial domain by starting from an established tumor-virus reaction-diffusion model, performing a center-manifold reduction that incorporates nonlinear terms to derive a complex Ginzburg-Landau amplitude equation, and estimating its parameters directly from the original kinetics. This reduced normal form equation explains the emergence of experimentally observed patterns—such as hollow rings and target waves—and shows that, at longer timescales, these patterns naturally evolve toward spiral waves and a turbulent regime. Our work provides a mechanistic link between the kinetic Hopf bifurcation and the rich spatiotemporal structures observed in oncolytic virotherapy models, suggesting that these patterns are not numerical artifacts but an intrinsic feature of the system.

1 Introduction

Cancer, defined by the uncontrolled growth and spread of abnormal cells, remains one of the leading causes of death worldwide [1]. Over the years, diverse treatments have been devised—chemotherapy [2], radiotherapy [3], immunotherapy [4], and more recently, oncolytic virotherapy [5]. Oncolytic virotherapy employs replicating viruses, which selectively infect and destroy tumor cells while sparing healthy tissues [6]. It looks very appealing in principle but in practice there are countless uncertainties at the biological level [7]. For example, how does the immune system react to this foreign virus [8]? How easily can tumor develop viral resistance [9]? These type of questions arise for any treatment in the oncologist's arsenal, but fortunately these uncertainties can be made more tangible by formulating them as mathematical questions, providing at least a benchmark for expected outcomes and make an educated guess about directions, which promise the most improvement [10].

If we adopt an ecological picture, we have oncolytic virus as a parasite trying to invade the ecosystem of the tumor host [11]. Host-parasite systems have been studied extensively, particularly in the models of epidemiological outbreaks [12] which highlight the basic reproduction number (R₀) as the main parameter to measure the contagiousness and the effect of the virus [13]. For oncolytic virotherapy, a similar approach has been taken to formulate Ordinary Differential Equations (ODEs) [14, 15] which identify a replication threshold beyond which the tumor population is expected to show oscillations. While such oscillations come naturally by a high enough replication rate, they have also been modeled through the lens of the viral lytic cycle [16, 17], where a time delay is used to model the recurrent viral bursts. In this paper however, we start from a previous model [10], which demonstrated that for high enough viral replication, these viral infection oscillations can drive the tumor to near-extinction. We note however that current viruses are still far from this Hopf bifurcation threshold [10]. Therefore, the natural question to ask here is: even if virus replication could be increased as much as desired, would that be enough for effective tumor eradication [18]?

An important aspect of the viral infection of cancer is the spatial extent of the tumor. Tumors have a finite size and the virus cannot reach all cells instantly, it needs time to infect, produce more virus, infect again etc. As a result periodic dynamics at different location along the tumor mix and exchange, leading to interesting spatio-temporal patterns [10]. Moreover, both in vitro and in silico approaches [9, 19] have demonstrated that the spatial structure of the tumor is crucial for therapeutic success. In particular, [19] used spatial models to show that the dynamics driven by waves of virus chasing the tumor lead to three distinct outcomes—hollow rings (successful clearance), concentric rings (tumor regrowth), and fragmented aggregates.

The previous model by Baabdulla et al. [10] consisted of three components - uninfected cancer, infected cancer and free virus (C, I, and V respectively) where virus diffuses more than the tumor due to its smaller size [10]. By simulating this Reaction-Diffusion system, the patterns observed before by [19] could be reproduced. In this work, we aim to demonstrate that these patterns are robust and not merely artifacts of some specific parameter choices. Our approach is to use mathematical reductions to make clear the basic mechanism underlying these patterns. We present a flow chart for our reduction steps in Figure 1.

Figure 1

Figure 1. Flow of analytical reductions from center-manifold dynamics to amplitude equations, parameter estimation, and pattern formation. The steps related to the well mixed ODE system are highlighted in yellow, while the results for the spatial PDE are shown with blue background.

In this approach, we first reduce the temporal dynamics of the full model to its center manifold [20] just above the Hopf bifurcation threshold. This reduction leads to a two-dimensional system, which can be related back to the original model and whose linear dynamics are directly given by the complex pair of eigenvalues. Further, to capture the resultant limit cycle in this reduced system, we determine the necessary nonlinear terms for the center manifold dynamics, which turn out to be quadratic since the original model's nonlinearities are also purely quadratic. However, since these nonlinear terms are coordinate dependent, we perform a normal form transformation through a lemma by Kuznetzsov [21] which allows us to interpret the reduced system in terms of the Stuart-Landau Equation (SLE) [22]. A similar framework involving normal form reductions has been recently applied to determine the stability of periodic solutions in delayed lytic cycle models [23]. This reduction constitutes the first analytical part of the paper, reiterating the well-known result that systems near a supercritical Hopf bifurcation can be universally represented by the SLE via normal form theory [24]. Upon non-dimensionalization, this results in a single asynchrony parameter β [25], which is a measure of frequency shifts between oscillations of varying amplitudes, and which we estimate numerically for our system.

After completely understanding the kinetic behavior in its normal form (the SLE), we proceed to the second part of this work, where the role of diffusion is concerned. This spatiotemporal extension gives the well-known Complex Ginzburg-Landau Equation (CGLE) [26] as normal form. The CGLE is known to generate a spectacular range of dynamic patterns, including plane waves, spiral waves, and spatiotemporal turbulence [25]. Spiral waves, for instance, are among the most striking patterns observed across diverse systems, from the Belousov-Zhabotinsky reaction [27] to electrical activity in cardiac tissue [28]. However, our numerical estimate of the parameter β suggests that the C–I–V system operates well beyond the regime that can support stable spiral waves. Consequently, turbulent patterns appear to be the norm rather than the exception for a large range of biological parameters in our system.

This work proceeds as follows: Section 2 introduces the C–I–V model and recalls the bifurcation structure as derived in the literature [10, 18, 29]. In Section 3 we focus on the ODEs of the well mixed version (no diffusion) and reduce the system near the Hopf bifurcation on the 2-dimensional center manifold to a complex Stuart-Landau equation. In Section 4 we add back diffusion and explain how the diffusion terms are carried through the reduction steps, leading to the complex Ginzburg-Landau equation (CGLE) as normal form. We show some numerical simulations for the full (C, I, V) model, followed by some simulations for spiral waves and their transition to turbulence for the CGLE. Finally, in Section 5, we conclude by discussing the biological implications of the turbulent behavior, the role for other factors such as immune system and resistance, and offer a heuristic explanation for periodic spot splitting phenomenon observed in an earlier work [10].

2 Model description and oscillations via Hopf bifurcation

Oncolytic viruses are engineered to exploit the differential susceptibility of tumor and normal cells. To capture the essential feedback between viral replication and tumor regrowth, we employ the nondimensionalized C–I–V model of [10], which extends the framework of [18]. The model describes three interacting densities: C(t, x): uninfected tumor cells, I(t, x): infected tumor cells and V(t, x): free virus particles. Their dynamics follow a system of reaction–diffusion equations:

$\begin{array}{ll}\frac{\partial C}{\partial t}=D_c\Delta C+C(1-C-I)-CV,& (\mathrm{1a})\end{array}$

$\begin{array}{ll}\frac{\partial I}{\partial t}=D_i\Delta I+CV-aI,& (\mathrm{1b})\end{array}$

$\begin{array}{ll}\frac{\partial V}{\partial t}=\Delta V+\theta I-\gamma V,& (\mathrm{1c})\end{array}$

where Δ denotes the Laplace operator, and D_c, D_i are the diffusion coefficients for cancer cells and infected cancer cells, respectively. Tumor proliferation is described through a logistic term C(1 − C − I) limited by space and resources. Viral infection is described by a mass-action term −CV, infected cells lyse at rate a, and produce new virions with rate θ. Viral particles are cleared with rate γ. Most parameters are normalized to equal 1, except the death rate of infected cells a, the removal rate of free virus particles γ, and the burst size θ. The dimensional version of the model, together with realistic parameter estimates and sensitivity analysis can be found in the cited literature [10, 18]. Together, these terms form a minimal yet complete model of virus infection-mediated tumor control, balancing growth, lysis, and diffusion.

It should be noted that we exclude two more terms, which are often included in these oncolytic virus models, see for example [29]. We ignore a viral loss term of the form −β₂CV in the last equation, since we assume, as in [10, 30–32], that the loss of virus particles due to an infection of a cell is minimal. Moreover, we do not include a natural death term for cancer cells. A death term in the C equation will shift the value of the Hopf bifurcation point, but it will not change the dynamics near that Hopf point.

We study Equation 1 on a two—dimensional domain Ω ⊂ ℝ² with zero—flux boundaries,

$\nabla C·\text{n}=\nabla I·\text{n}=\nabla V·\text{n}=0\text{\hspace{1em}on\hspace{1em}}\partial \Omega ,$

so that both cells and virions remain confined within Ω. This closed-domain setting mirrors in-vitro tumor spheroid experiments and petri-dish experiments [33], where viral spread occurs within a constrained microenvironment.

We follow [10] in adopting the values listed in Table 1. Among these, the viral production rate θ acts as the primary bifurcation parameter controlling infection strength and oscillatory onset.

Table 1

Table 1. Non-dimensional parameter set used in [10].

2.1 Kinetic behavior

To understand the mechanisms driving oscillations, we first suppress spatial dependence and analyze the corresponding well-mixed system:

$\begin{array}{ll}Ċ=C(1-C-I)-CV,& (\mathrm{2a})\end{array}$

$\begin{array}{ll}İ=CV-aI,& (\mathrm{2b})\end{array}$

$\begin{array}{ll}\dot{V}=\theta I-\gamma V.& (\mathrm{2c})\end{array}$

This reduction highlights the feedback loop between viral amplification and tumor recovery. System (Equation 2) has been analyzed previously, see for example in [10, 18] and we summarize the main properties here. The system admits three equilibria:

$\begin{array}{l}{E}_0=(0,0,0),E_1=(1,0,0),\end{array}$

$\begin{array}{ll}{E}_{+}=\left(C^{*},\frac{\gamma (1-C^{*})}{\theta +\gamma },\frac{\theta (1-C^{*})}{\theta +\gamma }\right),\text{\hspace{1em}where\hspace{1em}}{C}^{*}=\frac{a\gamma }{\theta }.& (3)\end{array}$

The condition for existence of the mixed equilibrium E₊ is C^* < 1 or θ > aγ, and E₊ state is born through a transcritical bifurcation at θ_T = aγ, as demonstrated in [10]. Biologically, this threshold marks the point where viral replication compensates for clearance, allowing infection to persist. As θ increases, E₊ eventually loses stability through a Hopf bifurcation at θ_H, producing sustained oscillations in tumor and viral populations. We show the bifurcation diagram in Figure 2, where θ is the bifurcation parameter and all other parameters are fixed as in Table 1. For three values of θ we show the time dynamics of C(t) in Figure 3.

Figure 2

Figure 2. Bifurcation diagram vs. θ for Equation 2: creation of the coexistence state at the transcritical threshold θ_T = 44.39 and loss of equilibrium stability at the Hopf threshold θ_H = 338.45.

Figure 3

Figure 3. Temporal dynamics of C(t) across parameter regimes of Equation 2. (a) tumor persistence (θ = 40); (b) stable coexistence (θ = 100); (c) settles down to periodic oscillations (θ = 360).

The oscillations observed in Figure 3 originate from a loss of stability of the coexistence equilibrium E₊. To determine the mechanism, we linearize the kinetic system (Equation 2) about E₊ and examine the eigenvalues of its Jacobian matrix:

$\begin{array}{ll}\text{J}(E_{+})=\left(\begin{array}{ccc}-C^{*}& -C^{*}& -C^{*}\\ \frac{\theta (1-C^{*})}{\theta +\gamma }& -a& C^{*}\\ 0& \theta & -\gamma \end{array}\right).& (4)\end{array}$

The characteristic polynomial of J(E₊) depends smoothly on the viral production rate θ. At a critical value θ = θ_H, one pair of complex conjugate eigenvalues λ_2,3 = ±iω_0H crosses the imaginary axis, while the remaining real eigenvalue λ₁ stays strictly negative. [10] proved that this crossing satisfies the standard conditions of a forward Hopf bifurcation.

Moreover, near θ_H, the system's behavior is governed by the two-dimensional eigenspace associated with the critical pair λ_2,3, while the remaining direction is rapidly damped (λ₁ < 0.). This separation of timescales allows reduction to a planar oscillatory subsystem on the corresponding center manifold.

3 Analysis and model reductions

We perform a systematic model reduction to obtain a Stuart-Landau normal form. We proceed as follows:

3.1 Reduction to the center manifold

At the Hopf threshold, the Jacobian J(E₊(θ_H)) possesses a purely imaginary eigenpair ±iω₀ and a single stable eigenvalue λ₁(θ_H) < 0. Let e₃ denote the stable eigenvector and $({\text{e}}_1,{\text{e}}_2={\overline{\text{e}}}_1)$ the conjugate pair spanning the oscillatory subspace. The dynamics decompose accordingly: fast decay along e₃ and sustained oscillations within span{e₁, e₂}. To make this decomposition explicit, we first shift the equilibrium to the origin:

$\begin{array}{ll}\tilde{C}=C-C^{*},Ĩ=I-I^{*},Ṽ=V-V^{*}.& (5)\end{array}$

Perturbations $(\tilde{C},Ĩ,Ṽ)$ are then expressed in the eigenspace basis as span{e₁, e₂} to leading order as

$\begin{array}{ll}(\tilde{C},Ĩ,Ṽ)=z(t){\text{e}}_1+\overline{z}(t){\overline{\text{e}}}_1=2\text{Re}\left[z(t){\text{e}}_1\right],& (6)\end{array}$

introducing the complex amplitude z(t) = u(t) + iv(t). The real variables (u, v) represent the reduced coordinates on the center manifold and form the foundation of the Stuart–Landau reduction carried out later.

We make this transformation explicit, by computing the eigenvectors e₁ and e₂. We normalize the eigenvector e₁ so that its third component is real and equal to 1. Then it has the general form

$\begin{array}{ll}{\text{e}}_1={(c_1+i{c}_2,c_3+i{c}_4,1)}^{\top },c_j\in ℝ,& (7)\end{array}$

where the coefficients c₁, c₂, c₃, c₄ will be found later. Substituting Equation 7 into Equation 6 gives

$\begin{array}{ll}\tilde{C}=(u+iv)(c_1+i{c}_2)+(u-iv)(c_1-i{c}_2)=2(c_{1}u-c_{2}v),& (8)\end{array}$

$\begin{array}{ll}Ĩ=(u+iv)(c_3+i{c}_4)+(u-iv)(c_3-i{c}_4)=2(c_{3}u-c_{4}v),& (9)\end{array}$

$\begin{array}{ll}Ṽ=(u+iv)+(u-iv)=2u.& (10)\end{array}$

which defines the forward mapping $(u,v)\mapsto (\tilde{C},Ĩ,Ṽ)$ with

$\begin{array}{ll}\tilde{C}=2(c_{1}u-c_{2}v),& (\mathrm{11a})\end{array}$

$\begin{array}{ll}Ĩ=2(c_{3}u-c_{4}v),& (\mathrm{11b})\end{array}$

$\begin{array}{ll}Ṽ=2u,& (\mathrm{11c})\end{array}$

Solving for (u, v) yields

$\begin{array}{ll}u=\frac{Ṽ}{2},v=\frac{c_3\tilde{C}-c_1Ĩ}{2(c_1{c}_4-c_2{c}_3)}.& (12)\end{array}$

Biologically, u corresponds to the viral component, while v represents a weighted average of the tumor components.

To determine e₁ explicitly, we solve

$\begin{array}{ll}\left(i{\omega }_{0H}\text{I}-\text{J}(E_{+})\right){\text{e}}_1=0,& (13)\end{array}$

with J(E₊) given in Equation 4. Setting ${\text{e}}_1={(l_1,l_2,1)}^{\top }$ gives the linear system

$\begin{array}{ll}\left(\begin{array}{ccc}i{\omega }_{0H}+C^{*}& C^{*}& C^{*}\\ -\frac{{\theta }_H\left(1-C^{*}\right)}{{\theta }_H+\gamma }& i{\omega }_{0H}+a& -C^{*}\\ 0& -{\theta }_H& i{\omega }_{0H}+\gamma \end{array}\right)\left(\begin{array}{c}{l}_1\\ l_2\\ 1\end{array}\right)=0.& (14)\end{array}$

From the third row we obtain

$\begin{array}{ll}{l}_2=\frac{\gamma +i{\omega }_{0H}}{{\theta }_H},& (15)\end{array}$

and substituting into the first row gives

$\begin{array}{ll}{l}_1=-\frac{C^{*}}{{\theta }_H}\frac{\gamma +{\theta }_H+i{\omega }_{0H}}{C^{*}+i{\omega }_{0H}}.& (16)\end{array}$

Thus,

$\begin{array}{ll}{\text{e}}_1=\left(\begin{array}{c}-\frac{C^{*}}{{\theta }_H}\frac{\gamma +{\theta }_H+i{\omega }_{0H}}{C^{*}+i{\omega }_{0H}}\\ \frac{\gamma +i{\omega }_{0H}}{{\theta }_H}\\ 1\end{array}\right).& (17)\end{array}$

Writing l₁ = c₁ + ic₂ and l₂ = c₃ + ic₄ yields the real coefficients

$\begin{array}{lll}{c}_1=-\frac{C^{*}}{{\theta }_H}\frac{(\gamma +{\theta }_H)C^{*}+{\omega }_{0H}^2}{C^{*2}+{\omega }_{0H}^2},& c_2=\frac{C^{*}{\omega }_{0H}}{{\theta }_H}\frac{\gamma +{\theta }_H-C^{*}}{C^{*2}+{\omega }_{0H}^2}.& (\mathrm{18a})\end{array}$

$\begin{array}{lll}{c}_3=\frac{\gamma }{{\theta }_H},& c_4=\frac{{\omega }_{0H}}{{\theta }_H}.& (\mathrm{18b})\end{array}$

Substituting into Equation 12 provides the full coordinate transformation between the physical variables (C, I, V) and the reduced oscillatory variables (u, v). In particular, the denominator in the expression for v is

$\begin{array}{ll}2(c_1{c}_4-c_2{c}_3)=-\frac{2C^{*}{\omega }_{0H}}{{\theta }_H^2(C^{*2}+{\omega }_{0H}^2)}\left({\omega }_{0H}^2+(\gamma +C^{*}){\theta }_H+{\gamma }^2\right),& (19)\end{array}$

which is negative at the Hopf bifurcation. Since the quantity described in Equation 19 is the denominator involved in transformation given by Equation 12, it being non-zero confirms that the transformation is well defined in a neighborhood around θ = θ_H. We also note that c₁ < 0 and c₂, c₃, c₄ > 0.

3.2 Nonlinear dynamics on the center manifold

For θ near the Hopf threshold θ_H we write the eigenvalue associated with e₁ in the general form

$\begin{array}{ll}\mu (\theta )={\lambda }_0(\theta )+i{\omega }_0(\theta ),{\lambda }_0({\theta }_H)=0,{\omega }_0({\theta }_H)={\omega }_{0H}>0.& (20)\end{array}$

The real part λ₀(θ) governs the growth or decay of oscillations, while ω₀(θ) sets their frequency. On the center manifold, using the complex coordinate z = u + iv from Equation 6, the leading-order linear dynamics take the form

$\begin{array}{ll}ż=({\lambda }_0+i{\omega }_0)z+O(|z{|}^2).& (21)\end{array}$

We show that the nonlinear terms are purely quadratic by a standard application of center manifold theory [20].

Lemma 1. The dynamics on the center manifold at the Hopf point are given by the nonlinear complex amplitude equation

$\begin{array}{ll}ż=({\lambda }_0+i{\omega }_0)z+(\frac{g_{20}}{2}{z}^2+g_{11}z\bar{z}+\frac{g_{02}}{2}{\bar{z}}^2),& (22)\end{array}$

with coefficients as defined in Equations 26, 27, 31.

Proof. The original C–I–V model (Equation 2) contains only quadratic nonlinearities (terms such as C², CI, CV). After translation to E₊, the constant terms vanish, leaving a purely quadratic remainder. Since the change of variables to (u, v) in Equation 12 is linear, this degree is preserved.

Hence, the reduced vector field on the center manifold is quadratic:

$\begin{array}{ll}ż=({\lambda }_0+i{\omega }_0)z+({\mathcal{P}}_1(u,v)+i{\mathcal{P}}_2(u,v)),& (23)\end{array}$

with

$\begin{array}{ll}\mathcal{P}(u,v)=\left(\begin{array}{c}{p}_1{u}^2+p_{2}uv+p_3{v}^2\\ q_1{u}^2+q_{2}uv+q_3{v}^2\end{array}\right),& (24)\end{array}$

for some coefficients p_i, q_i ∈ ℝ.

Expanding u², v², and uv in powers of z = u + iv and $\bar{z}=u-iv$ gives

$\begin{array}{ll}{u}^2=\frac{1}{4}(z^2+2z\bar{z}+{\bar{z}}^2),& (\mathrm{25a})\end{array}$

$\begin{array}{ll}{v}^2=\frac{1}{4}(-z^2+2z\bar{z}-{\bar{z}}^2),& (\mathrm{25b})\end{array}$

$\begin{array}{ll}uv=\frac{1}{4}(-i{z}^2+i{\bar{z}}^2),& (\mathrm{25c})\end{array}$

which leads to the compact complex form (Equation 22). with coefficients

$\begin{array}{ll}{g}_{20}=\frac{1}{2}[(p_1-p_3)-i{p}_2+q_2+i(q_1-q_3)],& (\mathrm{26a})\end{array}$

$\begin{array}{ll}{g}_{11}=\frac{1}{2}[(p_1+p_3)+i(q_1+q_3)],& (\mathrm{26b})\end{array}$

$\begin{array}{ll}{g}_{02}=\frac{1}{2}[(p_1-p_3)+i{p}_2-q_2+i(q_1-q_3)].& (\mathrm{26c})\end{array}$

It remains to determine the coefficients p_i and q_i at the Hopf point. By construction, $u=\frac{1}{2}Ṽ$, and from Equation 2 the $\dot{V}$-equation is purely linear,

$\dot{V}=\theta I-\gamma V.$

Using Ṽ = 2u and Ĩ = 2(c₃u − c₄v) from Equation 11 gives

$\dot{u}=\frac{1}{2}(\theta Ĩ-\gamma Ṽ)=(\theta c_3-\gamma )u-\theta c_{4}v.$

Substituting c₃ = γ/θ and c₄ = ω_0H/θ yields $\dot{u}=-{\omega }_{0H}v$, confirming that

$\begin{array}{ll}{p}_1=p_2=p_3=0.& (27)\end{array}$

To find the coefficients q₁, q₂, q₃, we consider the parametrization of ${\text{e}}_1={(c_1+i{c}_2,c_3+i{c}_4,1)}^{\top }$. Recall from Equation 11 that

$\tilde{C}=2(c_{1}u-c_{2}v),Ĩ=2(c_{3}u-c_{4}v),Ṽ=2u,$

and from Equation 12,

$v=\frac{c_3\tilde{C}-c_1Ĩ}{2\delta },\delta :=c_1{c}_4-c_2{c}_3\ne 0.$

The quadratic parts of the ODE system Equation 2 in shifted variables are

$\begin{array}{ll}{\dot{\tilde{C}}}_{\text{quad}}=-({\tilde{C}}^2+\tilde{C}Ĩ+\tilde{C}Ṽ),{\dot{Ĩ}}_{\text{quad}}=\tilde{C}Ṽ,{\dot{Ṽ}}_{\text{quad}}=0.& (28)\end{array}$

Because c₁, c₃, δ are constants,

$\begin{array}{ll}{\dot{v}}_{\text{quad}}=\frac{c_3{\dot{\tilde{C}}}_{\text{quad}}-c_1{\dot{Ĩ}}_{\text{quad}}}{2\delta }.& (29)\end{array}$

Expressing ${\tilde{C}}^2$, $\tilde{C}Ĩ$, and $\tilde{C}Ṽ$ in terms of (u, v) and simplifying yields

$\begin{array}{l}{\dot{\tilde{C}}}_{\text{quad}}=-4[(c_1^2+c_1{c}_3+c_1)u^2+(-2c_1{c}_2-(c_1{c}_4+c_2{c}_3)\\ -c_2)uv+(c_2^2+c_2{c}_4)v^2],& (\mathrm{30a})\end{array}$

$\begin{array}{ll}{\dot{Ĩ}}_{\text{quad}}=4(c_1{u}^2-c_{2}uv).& (\mathrm{30b})\end{array}$

Substituting Equation 30 into Equation 29 and matching terms with ${\dot{v}}_{\text{quad}}=q_1{u}^2+q_{2}uv+q_3{v}^2$ yields

$\begin{array}{ll}\begin{array}{l}{q}_1=-\frac{2c_1(1+c_3)(c_1+c_3)}{c_1{c}_4-c_2{c}_3},\\ q_2=\frac{2(c_1{c}_2+2c_1{c}_2{c}_3+c_1{c}_3{c}_4+c_2{c}_3^2+c_2{c}_3)}{c_1{c}_4-c_2{c}_3},\\ q_3=-\frac{2c_2{c}_3(c_2+c_4)}{c_1{c}_4-c_2{c}_3}.\end{array}& (31)\end{array}$

Remark 1. The lemma stated above holds for all values of θ, since it follows from the quadratic nonlinearity of the original system. However we note that p_i and q_i depend on the parameter θ, and to find them for θ ≠ θ_H would be more involved.

3.3 Normal form reduction

To pass from the quadratic center-manifold system (Equation 22) to a universal amplitude description, we remove quadratic terms in Equation 22 by a parameter-dependent near-identity change of variables, following the Hopf normal-form construction (cf. [21], Lemma 3.4) with bifurcation parameter θ (throughout, the coefficients p_i(θ), q_i(θ) and hence g_ij(θ) depend smoothly on θ).

Lemma 2 (Transformation to the Hopf normal form). Consider the complex amplitude equation on the center manifold

$\begin{array}{ll}ż=\mu (\theta )z+\frac{g_{20}(\theta )}{2}{z}^2+g_{11}(\theta )z\bar{z}+\frac{g_{02}(\theta )}{2}{\bar{z}}^2,& (32)\end{array}$

where μ(θ) = λ₀(θ)+iω₀(θ) with λ₀(θ_H) = 0 and ω₀(θ_H) = ω_0H > 0. There exists an invertible, θ-dependent near-identity change of complex coordinate

$\begin{array}{ll}z=y+\frac{h_{20}}{2}{y}^2+h_{11}yȳ+\frac{h_{02}}{2}{ȳ}^2,& (33)\end{array}$

valid for sufficiently small |y| and |θ−θ_H|, which transforms Equation 32 into an equation without quadratic terms:

$\begin{array}{ll}ẏ=\mu (\theta )y+O(|y{|}^3).& (34)\end{array}$

Proof. To obtain a normal form of a Hopf bifurcation we remove the quadratic terms in Equation 32 by a near-identity change of variables (Equation 33), where the constants h₂₀, h₁₁, h₀₂ will be determined at a later stage of the proof.

We begin from the near-identity change of variables defined in Equation 33, which upon inverting formally gives

$\begin{array}{ll}y=z-\frac{h_{20}}{2}{y}^2-h_{11}yȳ-\frac{h_{02}}{2}{ȳ}^2.& (35)\end{array}$

Since z = y + O(|y|²), we may replace each quadratic monomial in y by its z–counterpart up to cubic error:

$\begin{array}{ll}{y}^2=z^2+O(|z{|}^3),& (36)\end{array}$

$\begin{array}{ll}yȳ=z\bar{z}+O(|z{|}^3),& (37)\end{array}$

$\begin{array}{ll}{ȳ}^2={\bar{z}}^2+O(|z{|}^3).& (38)\end{array}$

Substituting Equations 36–38 into Equation 35 yields

$\begin{array}{ll}y=z-\frac{h_{20}}{2}{z}^2-h_{11}z\bar{z}-\frac{h_{02}}{2}{\bar{z}}^2+O(|z{|}^3).& (39)\end{array}$

Next, we want to use Equation 32 to write dynamical equation for ẏ. We start by differentiating Equation 39:

$\begin{array}{ll}ẏ=ż-h_{20}zż-h_{11}(ż\bar{z}+z\dot{\bar{z}})-h_{02}\bar{z}\dot{\bar{z}}+O(|z{|}^3).& (40)\end{array}$

And next we substitute ż from Equation 32, keeping only terms up to quadratic order. Since

$ż=\mu z+O(|z{|}^2),\dot{\bar{z}}=\bar{\mu }\bar{z}+O(|z{|}^2),$

Equation 40 becomes

$\begin{array}{l}ẏ=\mu z+(\frac{g_{20}}{2}-\mu h_{20})z^2+(g_{11}-\mu h_{11}-\bar{\mu }{h}_{11})z\bar{z}\\ +(\frac{g_{02}}{2}-\bar{\mu }{h}_{02}){\bar{z}}^2+O(|z{|}^3).& (41)\end{array}$

Finally, using z = y + O(|y|²) and collecting quadratic terms in y, we obtain

$\begin{array}{l}ẏ=\mu y+\frac{1}{2}(g_{20}-\mu h_{20})y^2+(g_{11}-\bar{\mu }{h}_{11})yȳ\\ +\frac{1}{2}(g_{02}-(2\bar{\mu }-\mu )h_{02}){ȳ}^2+O(|y{|}^3).& (42)\end{array}$

Thus, by choosing

$\begin{array}{ll}{h}_{20}=\frac{g_{20}}{\mu },h_{11}=\frac{g_{11}}{\bar{\mu }},h_{02}=\frac{g_{02}}{2\bar{\mu }-\mu },& (43)\end{array}$

we cancel all quadratic terms and obtain the cubic normal form (Equation 34). At θ = θ_H we have μ(θ_H) = iω_0H ≠ 0 and $2\bar{\mu }({\theta }_H)-\mu ({\theta }_H)=-3i{\omega }_{0H}\ne 0$, so the denominators remain nonzero for |θ−θ_H| sufficiently small.

Consequently, applying Lemma 2–22 removes the quadratic terms and generates cubic terms resulting from the nonlinear transformation. Denoting the transformed coordinate by y, we get the following equation:

$\begin{array}{ll}ẏ=\mu (\theta )y+\frac{G_{30}}{6}{y}^3+\frac{G_{21}}{2}{y}^2ȳ+\frac{G_{12}}{2}y{ȳ}^2+\frac{G_{03}}{6}{ȳ}^3+O(|y{|}^4).& (44)\end{array}$

To move from this equation to the normal form description, we use another standard transformation method (c.f [21], Lemma 3.5) to derive the Stuart–Landau normal form as the correct approximation to center manifold. This is presented in Lemma 3, where we also proceed to analytically determine the Lyapunov coefficient, which demonstrates the supercritical nature of Hopf Bifurcation.

Lemma 3 (Normal form derivation). The cubic equation for y generated by Lemma 2,

$\begin{array}{ll}ẏ=\mu (\theta )y+{\displaystyle \sum _{j+k=3}}\frac{G_{jk}(\theta )}{j!k!}{y}^j{ȳ}^k+O(|y{|}^4),& (45)\end{array}$

can be transformed by an invertible, parameter-dependent change of complex coordinate

$\begin{array}{ll}y=w+\frac{H_{30}}{6}{w}^3+\frac{H_{21}}{2}{w}^2\bar{w}+\frac{H_{12}}{2}w{\bar{w}}^2+\frac{H_{03}}{6}{\bar{w}}^3,& (46)\end{array}$

for all sufficiently small |θ−θ_H|, into an equation with only one cubic term:

$\begin{array}{ll}ẇ=({\lambda }_0(\theta )+i{\omega }_0(\theta ))w+\kappa (\theta )w|w{|}^2+O(|w{|}^4).& (47)\end{array}$

Proof. The inverse transformation from y to w can be written by substituting w = y + O(|w|³) in Equation 46:

$\begin{array}{ll}w=y-\frac{H_{30}}{6}{y}^3-\frac{H_{21}}{2}{y}^2ȳ-\frac{H_{12}}{2}y{ȳ}^2-\frac{H_{03}}{6}{ȳ}^3+O(|y{|}^4).& (48)\end{array}$

Differentiating with respect to time and substituting ẏ from Equation 45:

$\begin{array}{l}ẇ=ẏ-\frac{H_{30}}{2}{y}^2ẏ-\frac{H_{21}}{2}(2yȳẏ+y^2\dot{ȳ})-\frac{H_{12}}{2}(ẏ{ȳ}^2+2yȳ\dot{ȳ})\\ -\frac{H_{03}}{2}{ȳ}^2\dot{ȳ}+O(|y{|}^4)\\ =\mu y+(\frac{G_{30}}{6}-\frac{\mu H_{30}}{2})y^3+(\frac{G_{21}}{2}-\mu H_{21}-\bar{\mu }\frac{H_{21}}{2})y^2ȳ\\ +(\frac{G_{12}}{2}-\frac{\mu H_{12}}{2}-\bar{\mu }{H}_{12})y{ȳ}^2+(\frac{G_{03}}{6}-\frac{\bar{\mu }{H}_{03}}{2}){ȳ}^3+O(|y{|}^4).& (49)\end{array}$

Substituting y in terms of w (from Equation 46) and collecting powers yields:

$\begin{array}{l}ẇ=\mu w+\frac{1}{6}(G_{30}-2\mu H_{30})w^3+\frac{1}{2}(G_{21}-(\mu +\bar{\mu })H_{21})w^2\bar{w}\\ +\frac{1}{2}(G_{12}-2\bar{\mu }{H}_{12})w{\bar{w}}^2+\frac{1}{6}(G_{03}-(3\bar{\mu }-\mu )H_{03}){\bar{w}}^3+O(|w{|}^4).& (50)\end{array}$

By setting H₃₀ = G₃₀/2μ, $H_{12}=G_{12}/2\bar{\mu }$, and $H_{03}=G_{03}/(3\bar{\mu }-\mu )$, we eliminate all cubic terms except $w^2\bar{w}$. The coefficient of the $w^2\bar{w}$ term is $\frac{1}{2}(G_{21}-(\mu +\bar{\mu })H_{21})$. This term could be removed if $H_{21}=G_{21}/(\mu +\bar{\mu })$, but since $\mu ({\theta }_H)+\bar{\mu }({\theta }_H)=0$ at the bifurcation point, the divisor vanishes. For this reason, $w^2\bar{w}$ term is called resonant term and to ensure a smooth transformation with respect to θ, we set H₂₁ = 0. So after all these substitutions in Equation 50, we obtain:

$\begin{array}{l}ẇ=\mu w+\frac{G_{21}}{2}{w}^2\bar{w}+O(|w{|}^4),\end{array}$

which is the Stuart–Landau form with Lyapunov coefficient κ(θ) = G₂₁/2.

We already know that this Hopf bifurcation is supercritical from the stable limit cycle observed numerically (Figures 2, 3) but can demonstrate this also by checking the sign of κ(θ_H). From another result from [21]:

$\begin{array}{l}\kappa ({\theta }_H)=\frac{i}{2{\omega }_{0H}}(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{1}{3}|g_{02}{|}^2)\to \text{Re}(\kappa ({\theta }_H))\\ =\frac{i}{2{\omega }_{0H}}\text{Im}(g_{20}{g}_{11}).\end{array}$

Recall Equation 26 by noting that p_i = 0 (Equation 27), so we can write:

$g_{20}=\frac{1}{2}(q_2+i(q_1-q_3))g_{11}=\frac{i}{2}(q_1+q_3),$

which gives:

$\begin{array}{ll}\text{Re}(\kappa ({\theta }_H))=\frac{i}{2{\omega }_{0H}}\frac{q_2}{2}\frac{i(q_1+q_3)}{2}=-\frac{1}{8{\omega }_{0H}}{q}_2(q_1+q_3).& (51)\end{array}$

Since we already know q_i in terms of eigenvector components c_i (Equation 31), this can be computed and we obtain q₂ < 0, q₁ + q₃ < 0 which gives negative value of κ confirming that the bifurcation is supercritical and a stable limit cycle emerges from it. We note that this normal form approach to determine bifurcation elements is quite general and specifically for a class of delayed oncolytic model, formalized in [23].

3.4 Phase and amplitude equations

We now reinterpret the Stuart–Landau normal form

$\begin{array}{ll}ẇ=({\lambda }_0+i{\omega }_0)w-(a_1+i{b}_1)|w{|}^{2}w+O(|w{|}^4),& (52)\end{array}$

where λ₀(θ_H) = 0 and ω₀(θ_H) = ω_0H > 0, and where c₁(θ) = −a₁ − ib₁ is the cubic normal form coefficient.

Writing the complex amplitude as w = re^iϕ with r ≥ 0 and separating real and imaginary parts yields

$\begin{array}{ll}ṙ=r({\lambda }_0-a_1{r}^2)+O(r^4),& (\mathrm{53a})\end{array}$

$\begin{array}{ll}\dot{\varphi }={\omega }_0-b_1{r}^2+O(r^3),& (\mathrm{53b})\end{array}$

which, to leading order, are the classical Λ–Ω equations [25, 34] with

$\Lambda (r)={\lambda }_0-a_1{r}^2,\Omega (r)={\omega }_0-b_1{r}^2.$

For a supercritical Hopf bifurcation (a₁ > 0), the amplitude equation admits a stable limit cycle of radius $r_{*}^2={\lambda }_0/a_1$, oscillating at frequency $\Omega (r_{*})={\omega }_0-b_1{r}_{*}^2$.

As θ crosses θ_H, λ₀ changes sign and r^* grows from 0 for $\theta ={\theta }_H^{+}$ which agrees with the original ODE description, where Hopf Bifurcation creates a stable limit cycle with increasing amplitude as we move a bit further from the bifurcation point. We consider only the supercritical case because for the subcritical case a₁ < 0 there should be unstable cycles, which are not observed in simulations.

To reformulate (Equation 53) in a simpler form, we first remove the trivial rotation at frequency ω₀, we introduce the co-rotating variable W(t):

$\begin{array}{ll}w(t)=W(t)e^{i{\omega }_{0}t}\to W(t)=w(t)e^{-i{\omega }_{0}t},& (54)\end{array}$

so that |w| = |W|. Differentiation and substitution into Equation 52 give

$\begin{array}{l}Ẇ+i{\omega }_{0}W=({\lambda }_0+i{\omega }_0)W-(a_1+i{b}_1)|W{|}^{2}W+O(|z{|}^4),\\ Ẇ={\lambda }_{0}W-(a_1+i{b}_1)|W{|}^{2}W+O(|W{|}^4).& (55)\end{array}$

Assuming a₁ > 0, we scale variables to normalize the linear and cubic terms by defining

$\begin{array}{ll}W=\alpha \rho ,\alpha =\sqrt{\frac{{\lambda }_0}{a_1}},& (\mathrm{56a})\end{array}$

$\begin{array}{ll}s={\lambda }_{0}t\Rightarrow \frac{d}{dt}={\lambda }_0\frac{d}{ds}.& (\mathrm{56b})\end{array}$

Using ${\alpha }^2={\lambda }_0/a_1$ and substituting into Equation 55 yields the normal form

$\begin{array}{ll}\frac{d\rho }{ds}=\rho -(1+i\beta )|\rho {|}^2\rho +O(|\rho {|}^4)& (57)\end{array}$

where

$\begin{array}{ll}\beta =\frac{b_1}{a_1}& (58)\end{array}$

captures the strength of amplitude–frequency coupling.

In this nondimensional form, the amplitude saturates at |ρ| = 1 while the imaginary coefficient β|ρ|² generates the amplitude-dependent frequency shift in periodic orbits near unstable point. β in particular serves as an asynchrony parameter in the reduced normal form, since a higher β implies that even slight variations in amplitude can change the frequency by a large amount. Effects of this will become apparent in the spatial context, where such sensitive oscillators when coupled by diffusion would fail to synchronize and β would emerge as the determining factor in stability of different possible patterns.

3.5 Evaluating the normal-form parameter β

The parameter β in Equation 57 quantifies the nonlinear frequency shift per unit amplitude growth and links the reduced dynamics to the underlying C–I–V kinetics. To express it in measurable terms, note that for θ > θ_H the stable limit cycle with radius r_* satisfies

$\begin{array}{ll}\Lambda (r_{*})=0\Rightarrow a_1=\frac{{\lambda }_0}{r_{*}^2}.& (59)\end{array}$

Its oscillation frequency,

$\begin{array}{ll}{\Omega }_{\text{osc}}(\theta )={\omega }_0-b_1{r}_{*}^2,b_1=\frac{{\omega }_0-{\Omega }_{\text{osc}}}{r_{*}^2},& (60)\end{array}$

then gives

$\begin{array}{ll}\beta (\theta )=\frac{b_1}{a_1}=\frac{{\omega }_0(\theta )-{\Omega }_{\text{osc}}(\theta )}{{\lambda }_0(\theta )}=\frac{{\omega }_0(\theta )-2\pi /T_{\text{osc}}(\theta )}{{\lambda }_0(\theta )}.& (61)\end{array}$

Thus β can be obtained directly by combining the linear eigenvalues (λ₀, ω₀) of the coexistence state with the nonlinear oscillation period T_osc(θ) of the full system.

This observation leads to the following workflow: For each θ near the Hopf threshold:

1. Compute the leading complex eigenvalue μ(θ) = λ₀(θ)+iω₀(θ) of the Jacobian (Equation 4).

2. Integrate the ODE (Equation 2) until transients decay and measure the oscillation period T_osc(θ) from successive peaks of C(t).

3. Evaluate β(θ) via Equation 61.

Example: For our base parameter set from Table 1 we show in Figure 4 the eigenvalue pair of E₊ as θ increases: λ₀(θ) crosses zero at θ_H, while ω₀(θ) varies slowly. Figure 5 illustrates the resulting oscillations of C(t) for representative values of θ above the threshold. Combining linear and nonlinear data gives the nearly constant β(θ) values in Table 2. The ratio ω₀ − Ω_osc grows roughly proportionally to λ₀, keeping β of order 40 with only weak dependence on θ. This suggests that the frequency-amplitude coupling is an intrinsic property of the underlying C–I–V kinetics rather than a parameter-specific artifact.

Figure 4

Figure 4. Leading eigenvalue components of the coexistence equilibrium E₊ as functions of viral production rate θ. The real part λ₀(θ) crosses zero at θ = θ_H, marking the Hopf bifurcation, while the imaginary part ω₀(θ) varies only weakly.

Figure 5

Figure 5. Limit cycle oscillations of tumor cell density C(t) above the Hopf bifurcation. Each panel shows the last 100 time units of a simulation with total duration T = 5, 000. Red markers indicate peaks used to compute T_osc and Ω_osc.

Table 2

Table 2. Comparison of Jacobian eigenvalues and measured oscillation properties for θ = 340–360.

4 Inclusion of diffusion

To add diffusion back into the above transformation and scaling arguments we recall the linear change of variables from Equation 12:

$\begin{array}{ll}u=\frac{Ṽ}{2},v=\frac{c_3\tilde{C}-c_1Ĩ}{2(c_1{c}_4-c_2{c}_3)},& (62)\end{array}$

where the coefficients c₁, c₂, c₃, c₄ are given in Equation 18. Thus, u depends primarily on the viral perturbation Ṽ, while v is a fixed linear combination of the cellular perturbations $\tilde{C}$ and Ĩ.

Since the viral diffusivity was normalized to unity (D_V = 1), u inherits diffusion coefficient 1. Similarly, v is a linear combination of the two cancer populations, which have a common diffusion coefficient D_c = D_i. Hence v inherits the same effective diffusivity D_c. In the (u, v) coordinates the diffusion matrix is diagonal:

$\begin{array}{ll}D=\left(\begin{array}{cc}1& 0\\ 0& D_c\end{array}\right).& (63)\end{array}$

Combining spatial coupling with the linear kinetics, we obtain on the Center Manifold that

$\begin{array}{ll}\frac{\partial }{\partial t}\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{cc}{\lambda }_0& -{\omega }_0\\ {\omega }_0& {\lambda }_0\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)+D{\nabla }^2\left(\begin{array}{c}u\\ v\end{array}\right)+\text{nonlinear terms}.& (\mathrm{64a})\end{array}$

This reduced reaction-diffusion system captures the linear oscillatory modes and their spatial coupling. The next step is to include the weak nonlinearities that saturate these oscillations and ultimately give rise to pattern formation.

4.1 Observed patterns and explanation

In the previous section we reduced the oncolytic virotherapy model to an oscillatory form, highlighting its proximity to a Hopf bifurcation. Once such a reduction is available, the system can be interpreted within the general theory of coupled oscillators, where diffusion and nonlinear interactions organize local oscillations into coherent spatial structures such as plane waves, spiral waves, and turbulent patterns [25, 28, 35].

Here, we illustrate the behavior of the original C − −I − −V model (Equation 1) using VisualPDE [36] simulations. We encourage readers to explore the interactive models via the links included in the figure captions. We explore the dynamics for various values of θ and describe the main phenomena that emerge. Throughout these experiments, we impose Neumann boundary conditions. Initial conditions consist of a homogeneous tumor mass (C, I, V) = (1, 0, 0) with localized viral inoculation spots, where the concentration is set to (C, I, V) = (1, 0, 10). Depending on the experiment, the computational 2-dimensional domain is chosen to be either square or circular with Neumann boundary conditions, as indicated in the figures. These numerical experiments serve as a reference point for the more systematic analysis to follow, where we return to the analytically reduced system to isolate and interpret specific features.

Initial Ring Formation and Wave Propagation. To observe the initial behavior upon inoculation, we take circular domain with radius 250 units, and θ = 350, with two spots of virus injected. Upon inoculation, the infection spreads in waves of virus chasing the tumor to form expanding circular rings that radiate outward from the source. As shown in Figure 6, these waves originate at the inoculation site, which for θ > θ_H, act as a pacemaker sustaining periodic emission of fronts. But over time, we see that these fronts collide, merge, and distort–marking the transition from ordered wave propagation to the irregular patterns to come in later stages. We observe some symmetry in the patterns of Figure 6, the cause of which seems that since we have inoculated the spots at the same time, they have same initial behavior and the disturbances created by them add up linearly so that patterns are expected to be symmetric for initial times along the line which goes between them.

Figure 6

Figure 6. Early dynamics of system (Equation 1) following localized viral inoculation (VisualPDE Simulation [36]). Domain radius = 250, θ = 350, other parameters from Table 1. We have plotted the tumor concentration, where gray regions are close to zero value while the red regions cover high tumor densities. (a) Viral introduction at a point source. (b) Expanding circular fronts form around the core. (c) As fronts widen and interact, curvature decreases and the system progresses toward complex, disordered motion. The simulation can be accessed here: https://visualpde.com/sim/?mini=lfAVqB3Z.

Spatiotemporal patterns beyond Hopf bifurcation To see the long term behavior of the system (Equation 1), we take a square domain of size 500 units and simulate with viral production rate θ = 500, which is well above the Hopf threshold (recall that θ_H ≈ 338). A different color scheme is adopted to make the contrast in tumor density more apparent, where blue region is the lowest, red is intermediate, orange is high and yellow is the highest. Representative snapshots in Figure 7 show that the domain is tiled by fragmented tumor aggregates whose collisions and mergers sustain long-term turbulent activity rather than relaxing to a simple fixed pattern.

Figure 7

Figure 7. Persistent spatiotemporal activity beyond Hopf bifurcation (θ = 500). Snapshots of the field show turbulent organization of tumor mass on a square domain with side 500 units (VisualPDE; [36]). The simulation can be accessed here: https://visualpde.com/sim/?mini=LjbG7NBb. (a) t = 500, (b) t = 750, (c) t = 1, 000.

Spiral fragments To identify some organizing structure within this turbulent behavior, we consider a close-up sequence in Figure 8, again for circular domain of radius 250 units and θ = 350. In this zoomed in snapshot at the boundary of the circular domain, we observe a single spiraling disturbance, which seems to rotate and also drift a little, with possibly a length scale for its wavelength and curvature. This serves as a visual clue to look for spiral wave solutions for our original system (Equation 1), which we will discuss in more detail in the upcoming sections.

Figure 8

Figure 8. Plot of tumor mass density showing closeup on spiral disturbance for θ = 350, with color intensity proportional to tumor concentration. Consecutive frames show a single rotating wave segment where the tip precesses slowly while the arm maintains a rough structure and length scale, characteristic of spiral-wave dynamics (adapted from VisualPDE [36]). The simulation can be accessed here: https://visualpde.com/sim/?mini=oXulWk5A.

Amplitude modulations at fixed spatial points To further probe the spatiotemporal behavior, for our choice of parameter θ = 350, we restrict attention to the temporal evolution at a fixed spatial point in the domain. We would like to choose a point not too close to the boundary, nor at the center, nor at a symmetrical point so we choose x₀ = (10, 20) for our purpose here and record the time series of V at x₀ for a simulation in square domain of side L = 100. The resulting dynamics, shown in Figure 9, reveal a slow envelope riding on top of the fast oscillations near the natural frequency. To better visualize the amplitude modulations, we need to tune out the fast oscillations for which we construct a Poincaré sampling at fixed intervals Δt close to one oscillation period (here for Δt = 4.2s). The resulting trajectory (Figure 9b) shows roughly how much the amplitude varies over successive return to the same phase of oscillation.

Figure 9

Figure 9. Local temporal structure at a fixed spatial point x₀ = (10, 20) for θ = 350 and domain size 100 × 100. On the x-axis, we have time and y-axis is V = V(10, 20, t). (a) A slow envelope modulates the fast oscillation near the natural frequency. (b) Sampled once per oscillation, points follow a thin curve clearly showing the amplitude modulations.

Phase portraits across θ Finally, instead of monitoring the temporal evolution at a single component, we consider the three-dimensional trajectory (C(x₀, t), I(x₀, t), V(x₀, t)) again with fixed x₀ = (10, 20). To get a clean picture which highlights amplitude modulation, we sample again near the natural frequency, as done for previous Figure 9 with Δt = 4.2s, to obtain the phase portraits shown in Figure 10. Note that here the axes are normalized so that the unstable fixed point E₊ sits at (1, 1, 1). Closer to the onset of Hopf bifurcation (θ = 344), the trajectory forms a relatively thin closed loop in (C, I, V); by θ = 350, the loop thickens into a torus-like annulus and the instantaneous amplitudes vary more widely.

Figure 10

Figure 10. Poincaré phase portraits in (C, I, V) obtained by sampling near the natural frequency at Δt = 4.2s. The x, y, z axes are respectively $(C,I,V)=(\frac{C(10,20,t)}{C^{*}},\frac{I(10,20,t)}{I^{*}},\frac{V(10,20,t)}{V^{*}})$, where we normalize the coordinates so that the unstable point E₊ sits at (1, 1, 1). Just beyond onset (θ_H ≈ 338), trajectories lie on a thin closed loop (approximate limit cycle). At higher value (θ = 350), the system shows signs of period doubling where the original limit cycle can shrink to smaller amplitude cycles and then make a return back to the limit cycle. (a) θ = 344 ≳ θ_H. (b) θ = 350 ≫ θ_H.

4.2 Explanation

The natural starting point for analysis is therefore to find a family of plane-wave solutions of the cubic complex Ginzburg–Landau equation (CGLE) [26] with real diffusion in two spatial dimensions:

$\begin{array}{ll}{\partial }_t\rho =\rho -(1+i\beta )|\rho {|}^2\rho +D\Delta \rho ,x\in {ℝ}^2,D>0,\beta \in ℝ,& (65)\end{array}$

where ρ(x, t) is the complex amplitude field over the spatial domain and Δ = ∂_xx + ∂_yy is the Laplacian.

However, we note that in matrix notation we had D_u, D_v for diffusion coefficients but for the amplitude equation, we have a single D. Here we take D = (D_u + D_v)/2 ≈ 0.5 which is supported by the work of [35], who showed that amplitude equations for plane waves and spiral waves are governed by an effective diffusion rate, which is the average of the two interacting components. But since our goal here is to just explain the patterns and not to make exact predictions, we will be more relaxed with the exact value of D such as in spiral wave simulations for next section, where we take D = 1.

The complex Ginzburg Landau equation shows a number of spatial patterns, which we consider in sequence.

4.2.1 Spatially uniform oscillations

The spatially uniform oscillation is represented as ρ(x, t) = e^−iβt, which is synchronized oscillation of the whole domain with unit amplitude and frequency β.

4.2.2 Plane waves

Take the plane-wave ansatz with amplitude α ∈ (0, 1] and wavevector k ∈ ℝ²,

$\begin{array}{ll}\rho (\text{x},t)=\alpha e^{i(\text{k}·\text{x}-\Omega t)}.& (66)\end{array}$

Then |ρ|² = α², ∂_tρ = −iΩρ, and Δρ = −|k|²ρ. Substituting into Equation 65 and dividing by ρ gives

$\begin{array}{ll}-i\Omega =1-(1+i\beta ){\alpha }^2-D|\text{k}{|}^2=(1-{\alpha }^2-D|\text{k}{|}^2)-i\beta {\alpha }^2.& (67)\end{array}$

Equating real and imaginary parts, we find

$\begin{array}{ll}\Omega =\beta {\alpha }^2,1-{\alpha }^2-D|\text{k}{|}^2=0.& (68)\end{array}$

Thus the admissible plane-wave family is

$\begin{array}{ll}\rho (\text{x},t)=\alpha exp(i\left[\text{k}·\text{x}-\beta {\alpha }^{2}t\right]),& (69)\end{array}$

for any wave vector k with

$|\text{k}|=\sqrt{\frac{1-{\alpha }^2}{D}}.$

The homogeneous oscillation without spatial propagation corresponds to α = 1, for which k = 0 and we obtain ρ(x, t) = e^−iβt.

Eliminating α from Equation 68 gives the dispersion relation between frequency Ω and wave mode k

$\begin{array}{ll}\Omega (\text{k})=\beta (1-D|\text{k}{|}^2),|\text{k}|\le D^{-1/2}.& (70)\end{array}$

Phase and group velocities are colinear with k:

$\begin{array}{ll}{\text{v}}_p=\frac{\Omega }{|\text{k}|}\widehat{\text{k}}=\beta (\frac{1}{|\text{k}|}-D|\text{k}|)\widehat{\text{k}},{\text{v}}_g={\nabla }_{\text{k}}\Omega =-2\beta D\text{k},& (71)\end{array}$

so v_p and v_g always point in opposite directions when k ≠ 0. We see that the magnitudes of the phase velocity v_p and the group velocity v_g increase linearly in the coefficient β of the cubic term. An increase in the diffusion coefficient D reduces the phase speed and increases the group velocity and as D → |k|⁻² the phase becomes stationary.

The plane waves can be seen in the circular spread patterns in Figure 6, where for larger radius, the waves moves in normal direction as a plane wave.

4.2.3 Spiral waves and turbulence

Spiral waves present analytic challenges that their plane wave counterparts do not. As noted earlier, spiral dynamics originate at a central core, corresponding to the limit r → 0 in polar coordinates centered on the core. This introduces singular terms that must be handled carefully to avoid divergences at the origin, while at the same time ensuring smooth continuity with the far–field plane–wave behavior. [37] developed techniques to bridge this core–to–tail connection, providing the first demonstration that reaction–diffusion dynamics alone can sustain spiral waves without the need for additional terms. Building on this, [35] formulated a general theory of oscillatory media near Hopf bifurcation, reducing such systems to the λ–ω class. As Murray's classic text explains in detail [38], this reduction ultimately leads to the complex Ginzburg–Landau equation (CGLE) (Equation 65), which admits both plane–wave and spiral solutions [25].

Comprehensive reviews, most notably the survey by [26], have since documented the full breadth of spiral-wave phenomena, including their stability, the mechanisms of spiral breakup, and transitions to turbulence governed by the cubic CGLE. Returning to Equation 65, the single control parameter β quantifies the strength of amplitude—frequency coupling. The types of spirals that can arise is huge, including linear and logarithmic spirals as well as one-armed, two- armed or m-armed spirals, and meandering sprials. All these are well documented in the cited literature [26, 35, 37–39] and we will not recall the underlying scaling arguments here. Instead, we reproduce the main spiral regimes through direct numerical simulations and compare them with the theoretical expectations.

4.3 Numerical simulation of spiral waves

To illustrate the spiraling solutions of Equation 65, we perform direct numerical simulations on a two–dimensional square domain, with Neumann boundary conditions. The domain size is L = 100 with N = 256 grid points in each direction, giving a spatial resolution Δx = L/N ≈ 0.39. The Laplacian is discretized using second–order central differences, and time integration is carried out with an explicit Euler scheme with timestep Δt = 10⁻³. To visualize the spiral waves we have chosen an explicit solver instead of VisualPDE, since an explicit solver allows us full control on the discretization, time stepping, and initial condition. Our choices above are sufficient to resolve the spiral dynamics without visible numerical artifacts.

To seed a spiral, we prescribe

$\begin{array}{ll}\rho (x,y,0)=tanh(2r)exp(i\theta +\epsilon (x,y)),& (72)\end{array}$

where $r=\sqrt{{(x-L/2)}^2+{(y-L/2)}^2}$ is the distance from the domain center and θ = arctan2(y − L/2, x − L/2) the polar angle. The prefactor tanh(2r) creates a smooth core, while the phase factor e^iθ winds the oscillation once around the center, imposing a single spiral seed. A small random perturbation $\epsilon (x,y)~\mathcal{U}(-10^{-3},10^{-3})$ is added to break perfect rotational symmetry and prevent numerical artifacts.

To interpret the spiral dynamics of Equation 65, we monitor both the amplitude |ρ(x, y, t)| and the phase arg(ρ(x, y, t)). In practice, we display the amplitude defect 1 − |ρ|² as a color map, which highlights spiral cores as bright spots, and overlay contours of constant phase to reveal the rotating arms and their interactions. This joint representation makes it possible to simultaneously follow the wavefront geometry and the location of topological defects. An example is shown in Figure 11.

Figure 11

Figure 11. Spiral wave dynamics for CGLE (Equation 65) for β = 1.5 and D = 1 at t = 80s in square domain of size 100 obtained from the seeded initial condition (Equation 72). Color shading shows the amplitude defect (1 − |ρ|²), making spiral cores visible as bright spots, while overlaid contours trace equal phase lines of the oscillation, revealing the structure of the rotating arms. The initial seed had one spiral at the center, but many more spiral centers arise over time.

4.4 Onset of turbulence

The dynamics of Equation 65 depend strongly on the amplitude–frequency coupling parameter β. For small values of |β|, the medium supports frozen spirals: coherent, steadily rotating structures with smooth arms (Figure 12a). Far from the spiral cor, the wavefronts match the plane–wave dispersion relation (Equation 70) derived in the previous section, while the core itself serves as a continuous source of rotational disturbances.

Figure 12

Figure 12. Spiral wave states in the CGLE (Equation 65) with D = 1. (a) At small β, the system sustains frozen spirals: stable rotating waves with coherent arms. (b) At larger β, nonlinear frequency shifts destabilize the arms, leading to spiral breakup and a turbulent regime with proliferating defects.

As |β| increases, nonlinear frequency shifts destabilize the spiral arms and it becomes harder to sustain synchronous oscillations across domains of differing amplitude. The wavefronts begin to develop irregularities, tips break apart, and topological defects proliferate, ultimately leading to defect turbulence (Figure 12b).

Here we refer to a more recent work by [39], which studied spiral-wave stability in spatial rock–paper–scissors models (which also reduce to the same CGLE we consider here). They identified three regimes: (i) the bound state (BS) phase of long-lived frozen spirals (small β, Figure 12a); (ii) the Eckhaus instability (EI) phase, where far-field perturbations destabilize spiral arms but the core remains intact; and (iii) the absolute instability (AI) phase, where instability propagates into the core and prevents any sustained spirals. Our simulations (Figure 12) reproduce this sequence—from frozen spirals at low β to turbulence at larger β. However, we cannot distinguish between the EI and AI phases in our numerical experiments.

5 Conclusion

To summarize, we started with a basic Reaction-Diffusion (RD) model for the replication and spread of an oncolytic virus in a tumor population, treating the viral burst size as the primary control parameter. We specifically wished to explore and explain the patterns that arise right above the Hopf bifurcation threshold, which we achieved through systematic reduction to the Complex Ginzburg-Landau Equation (CGLE) (Equation 65). In light of the normal form approach taken here, we identified the parameter β, which measures the asynchrony of oscillations (the degree to which frequency depends on amplitude), as the critical factor. It is known through previous work on amplitude-dependent frequency oscillators, and we also saw through our simulations, that this β term is responsible for the stability of different patterns in the medium.

Pattern formation in oncolytic virotherapy as we studied has been observed before in much more complicated models involving immune interactions [40], evolving resistance [9], agent-based simulations [19], detailed tumor geometries [41], and nuanced movement strategies (like pressure-driven flow) [42]. We note that if cellular movement is modeled by cross-diffusion, we can get not just spiral turbulence but also Turing instabilities [43] and . However, in this work we take a broader approach where for a treatment that is inherently oscillating, the implications on the spatiotemporal patterns can be identified by β, a parameter that is both simple to describe and easy to calculate. Specifically, if the eigenvalue at the unstable point is λ₀ + iω₀ and the frequency of the resulting limit cycle is ω, then β is proportional to the difference in frequencies relative to the distance from the threshold. The parameter β is proportional to $\frac{\omega -{\omega }_0}{{\lambda }_0}$. It is known that if β is smaller than a certain critical value (around 3, depending on the system), we can expect stable spiral patterns, while higher values of β lead to turbulent patterns [26] as the general outcome. However, [19] found a very curious result in their experiments, where the same system with the same initial conditions could lead to different treatment outcomes. They attributed this unpredictability either to intrinsic stochasticity within the model or to external anti-viral factors not included in the model [19], and the same problem persists in our model too, but it is also the natural side effect of reducing the dynamics too much.

Finally we reconsider the mathematical approach taken here, which reduced the complex dynamics to the picture of coupled nonlinear oscillators. This could be used to explain the main patterns already, but we can go further. In particular, a previous study [10] found that an inoculated spot of virus can split into two at regular intervals, a phenomenon termed periodic spot-splitting. We can now offer a possible explanation for this using the idea of amplitude-dependent frequency of the oscillators. In the initial inoculation, the center of the spot has the highest viral load (amplitude), while a point near the edge has a lower value. If the oscillators were truly linear, we would expect the entire spot to recover simultaneously after one oscillation. However, because we have shown the system is governed by the nonlinear Stuart-Landau Equation, the frequency of the highly infected spot center will be different from the less infected off-center points. This leads to a growing phase gap between the center and the periphery. This phase difference could create a scenario where the center of the spot reaches its minimum biomass (near-extinction) while some off-center point is still near its maximum–a desynchronization that would manifest visually as the splitting of the original inoculated spot. While this explanation remains heuristic, it clearly motivates a direction for further mathematical inquiry: establishing the precise conditions on the parameter β that determine the onset, choice of direction for symmetry breaking and the subsequent dynamics of spot splitting, thereby contributing to understanding of this new type of pattern formation for its own sake.

To conclude, we provide here a robust, minimal mathematical framework for understanding the failure modes of oncolytic virotherapy in spatially extended tumors. Our analysis demonstrates that, while constant efforts are made to increase viral contagiousness, high replication rates are mathematically destined to lead to fragmented patterns. This inherent instability ensures unreliable and incomplete viral coverage, confirming that maximizing viral replication still leads to therapeutic failure [18]. This finding supports the existing consensus that oncolytic viruses due to their dynamic limitations, cannot reliably clear a tumor entirely on their own and should therefore be strategically used in combination with other therapies [44–46]. This necessitates a shift from monotherapy optimization to designing viruses from scratch for synergistic action within combination treatment strategies.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://visualpde.com/.

Author contributions

TB: Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. TH: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. TH was supported through a Discovery Grant of the Natural Science and Engineering Research Council of Canada (NSERC), RGPIN-2023-04269.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | Crossref Full Text | Google Scholar

2. DeVita Jr VT, Chu E, A. history of cancer chemotherapy. Cancer Res. (2008) 68:8643–53. doi: 10.1158/0008-5472.CAN-07-6611

Crossref Full Text | Google Scholar

3. Gianfaldoni S, Gianfaldoni R, Wollina U, Lotti J, Tchernev G, Lotti T. An overview on radiotherapy: from its history to its current applications in dermatology. Open Access Maced J Med Sci. (2017) 5:521. doi: 10.3889/oamjms.2017.122

PubMed Abstract | Crossref Full Text | Google Scholar

4. Mellman I, Coukos G, Dranoff G. Cancer immunotherapy comes of age. Nature. (2011) 480:480–9. doi: 10.1038/nature10673

PubMed Abstract | Crossref Full Text | Google Scholar

5. Russell SJ, Peng KW, Bell JC. Oncolytic virotherapy. Nat Biotechnol. (2012) 30:658–70. doi: 10.1038/nbt.2287

PubMed Abstract | Crossref Full Text | Google Scholar

6. Lin D, Shen Y, Liang T. Oncolytic virotherapy: basic principles, recent advances and future directions. Signal Transd Target Therapy. (2023) 8:156. doi: 10.1038/s41392-023-01407-6

PubMed Abstract | Crossref Full Text | Google Scholar

7. Volovat SR, Scripcariu DV, Vasilache IA, Stolniceanu CR, Volovat C, Augustin IG, et al. Oncolytic virotherapy: a new paradigm in cancer immunotherapy. Int J Mol Sci. (2024) 25:1180. doi: 10.3390/ijms25021180

PubMed Abstract | Crossref Full Text | Google Scholar

8. Storey KM, Jackson TL. An agent-based model of combination oncolytic viral therapy and anti-PD-1 immunotherapy reveals the importance of spatial location when treating glioblastoma. Cancers. (2021) 13:5314. doi: 10.3390/cancers13215314

PubMed Abstract | Crossref Full Text | Google Scholar

9. Bhatt DK, Chammas R, Daemen T. Resistance mechanisms influencing oncolytic virotherapy, a systematic analysis. Vaccines. (2021) 9:1166. doi: 10.3390/vaccines9101166

PubMed Abstract | Crossref Full Text | Google Scholar

10. Baabdulla AA, Hillen T. Oscillations in a spatial oncolytic virus model. Bull Math Biol. (2024) 86:93. doi: 10.1007/s11538-024-01322-z

PubMed Abstract | Crossref Full Text | Google Scholar

11. Schoeps B, Lauer UM, Elbers K. Deciphering permissivity of human tumor ecosystems to oncolytic viruses. Oncogene. (2025) 44:1069–77. doi: 10.1038/s41388-025-03357-5

PubMed Abstract | Crossref Full Text | Google Scholar

12. Ezanno P, Vergu E, Langlais M, Gilot-Fromont E. Modelling the dynamics of host-parasite interactions: basic principles. In: New Frontiers of Molecular Epidemiology of Infectious Diseases. Cham: Springer (2011). p. 79–101.

Google Scholar

13. Delamater PL, Street EJ, Leslie TF, Yang YT, Jacobsen KH. Complexity of the basic reproduction number (R0). Emerg Infect Dis. (2019) 25:1. doi: 10.3201/eid2501.171901

PubMed Abstract | Crossref Full Text | Google Scholar

14. Bajzer Ž, Carr T, Josić K, Russell SJ, Dingli D. Modeling of cancer virotherapy with recombinant measles viruses. J Theoret Biol. (2008) 252:109–22. doi: 10.1016/j.jtbi.2008.01.016

PubMed Abstract | Crossref Full Text | Google Scholar

15. Tian JP. The replicability of oncolytic virus: defining conditions in tumor virotherapy. Mathem Biosci Eng. (2011) 8:841–60. doi: 10.3934/mbe.2011.8.841

PubMed Abstract | Crossref Full Text | Google Scholar

16. Najm F, Yafia R, Aziz-Alaoui M. Hopf bifurcation in oncolytic therapeutic modeling: Viruses as anti-tumor means with viral lytic cycle. Int J Bifurcat Chaos. (2022) 32:2250171. doi: 10.1142/S0218127422501711

Crossref Full Text | Google Scholar

17. Kayan Ş, Merdan H, Yafia R, Goktepe S. Bifurcation analysis of a modified tumor-immune system interaction model involving time delay. Mathem Model Nat Phenomena. (2017) 12:120–45. doi: 10.1051/mmnp/201712508

Crossref Full Text | Google Scholar

18. Pooladvand P, Yun CO, Yoon AR, Kim PS, Frascoli F. The role of viral infectivity in oncolytic virotherapy outcomes: a mathematical study. Math Biosci. (2021) 334:108520. doi: 10.1016/j.mbs.2020.108520

PubMed Abstract | Crossref Full Text | Google Scholar

19. Wodarz D, Hofacre A, Lau JW, Sun Z, Fan H, Komarova NL. Complex spatial dynamics of oncolytic viruses in vitro: mathematical and experimental approaches. PLoS Comput Biol. (2012) 8:e1002547. doi: 10.1371/journal.pcbi.1002547

PubMed Abstract | Crossref Full Text | Google Scholar

20. Carr J. Applications of Centre Manifold Theory, vol 35. Cham: Springer Science & Business Media. (2012).

Google Scholar

21. Kuznetsov YA. Elements of Applied Bifurcation Theory. Cham: Springer (1998).

Google Scholar

22. Landau LD. On the problem of turbulence. Dokl. Akad. Nauk USSR. (1944) 44:311.

Google Scholar

23. Najm F, Ahmed M, Yafia R, Aziz Alaoui M, Boukrim L. Hopf bifurcation and normal form in a delayed oncolytic model. Int J Biomathem. (2025) 18:2350111. doi: 10.1142/S1793524523501115

Crossref Full Text | Google Scholar

24. Guckenheimer J, Holmes P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, vol 42. Cham: Springer Science & Business Media. (2013).

Google Scholar

25. Kuramoto Y. Chemical Oscillations, Waves, and Turbulence. Cham: Springer Series in Synergetics. (1984). doi: 10.1007/978-3-642-69689-3

Crossref Full Text | Google Scholar

26. Aranson IS, Kramer L. The world of the complex Ginzburg-Landau equation. Rev Mod Phys. (2002) 74:99. doi: 10.1103/RevModPhys.74.99

Crossref Full Text | Google Scholar

27. Zaikin AN, Zhabotinsky AM. Concentration wave propagation in two-dimensional liquid-phase self-oscillating system. Nature. (1970) 225:535–7. doi: 10.1038/225535b0

PubMed Abstract | Crossref Full Text | Google Scholar

28. Winfree AT. Spiral waves of chemical activity. Science. (1972) 175:634–6. doi: 10.1126/science.175.4022.634

PubMed Abstract | Crossref Full Text | Google Scholar

29. Al-Tuwairqi SM, Al-Johani NO, Simbawa EA. Modeling dynamics of cancer virotherapy with immune response. Adv Differ Equ. (2020) 2020:438. doi: 10.1186/s13662-020-02893-6

Crossref Full Text | Google Scholar

30. Wodarz D. Viruses as antitumor weapons: defining conditions for tumor remission. Cancer Res. (2001) 61:3501–7.

PubMed Abstract | Google Scholar

31. Eftimie R, Dushoff J, Bridle BW, Bramson JL, Earn DJD. Multi-stability and multi-instability phenomena in a mathematical model of tumor-immune-virus interactions. Bull Math Biol. (2011) 73:2932–61. doi: 10.1007/s11538-011-9653-5

PubMed Abstract | Crossref Full Text | Google Scholar

32. Wu JT, Byrne HM, Kirn DH, Wein LM. Modeling and analysis of a virus that replicates selectively in tumor cells. Bull Math Biol. (2001) 63:731–68. doi: 10.1006/bulm.2001.0245

PubMed Abstract | Crossref Full Text | Google Scholar

33. Baabdulla A, Cristi F, Shmulevitz M, Hillen T. Mathematical modelling of reoviruses in cancer cell cultures. PLoS ONE. (2025) 13:e0190892. doi: 10.1371/journal.pone.0318078

PubMed Abstract | Crossref Full Text | Google Scholar

34. Kopell N, Howard L. Plane wave solutions to reaction-diffusion equations. Stud Appl Mathem. (1973) 52:291–328. doi: 10.1002/sapm1973524291

Crossref Full Text | Google Scholar

35. Duffy MR, Britton NF, Murray JD. Spiral wave solutions of practical reaction-diffusion systems. SIAM J Appl Math. (1980) 39:8–13. doi: 10.1137/0139002

Crossref Full Text | Google Scholar

36. Walker BJ, Townsend AK, Chudasama AK, Krause AL. VisualPDE: rapid interactive simulations of partial differential equations. Bull Math Biol. (2023) 85:113. doi: 10.1007/s11538-023-01218-4

PubMed Abstract | Crossref Full Text | Google Scholar

37. Cohen DS, Neu JC, Rosenblat S. Rotating spiral wave solutions of reaction-diffusion equations. SIAM J Appl Math. (1978) 35:536–47. doi: 10.1137/0135045

Crossref Full Text | Google Scholar

38. Murray JD. Mathematical Biology: I. An Introduction, vol 17. Cham: Springer Science & Business Media (2007).

Google Scholar

39. Szczesny B, Mobilia M, Rucklidge AM. Characterization of spiraling patterns in spatial rock-paper-scissors games. Physical Review E. (2014) 90:032704. doi: 10.1103/PhysRevE.90.032704

PubMed Abstract | Crossref Full Text | Google Scholar

40. Storey KM, Lawler SE, Jackson TL. Modeling oncolytic viral therapy, immune checkpoint inhibition, and the complex dynamics of innate and adaptive immunity in glioblastoma treatment. Front Physiol. (2020) 11:151. doi: 10.3389/fphys.2020.00151

PubMed Abstract | Crossref Full Text | Google Scholar

41. Bhatt DK, Janzen T, Daemen T, Weissing FJ. Modelling the spatial dynamics of oncolytic virotherapy in the presence of virus-resistant tumour cells. PLOS Comput Biol. (2022) 18:e1010076. doi: 10.1371/journal.pcbi.1010076

PubMed Abstract | Crossref Full Text | Google Scholar

42. Morselli D, Delitala ME, Frascoli F. Agent-based and continuum models for spatial dynamics of infection by oncolytic viruses. Bull Math Biol. (2023) 85:92. doi: 10.1007/s11538-023-01192-x

PubMed Abstract | Crossref Full Text | Google Scholar

43. Najm F, Yafia R, Aziz Alaoui MA. Turing bifurcation induced by cross-diffusion and amplitude equation in oncolytic therapeutic model: viruses as anti-tumor means. Int J Bifurcat Chaos. (2023) 33:2350062. doi: 10.1142/S0218127423500621

Crossref Full Text | Google Scholar

44. Wang Y, Zhu M, Chi H, Liu Y, Yu G. The combination therapy of oncolytic virotherapy. Front Pharmacol. (2024) 15:1380313. doi: 10.3389/fphar.2024.1380313

Crossref Full Text | Google Scholar

45. Fu R, Qi R, Xiong H, Lei X, Jiang Y, He J, et al. Combination therapy with oncolytic virus and T cells or mRNA vaccine amplifies antitumor effects. Signal Transd Target Therapy. (2024) 9:118. doi: 10.1038/s41392-024-01824-1

PubMed Abstract | Crossref Full Text | Google Scholar

46. Mohammadnejad N, Hillen T. Qualitative optimization of oncolytic virotherapy and immune therapy combination treatments. bioRxiv. 20250915676375. (2025). Available online at: https://www.biorxiv.org/content/early/2025/09/16/2025.09.15.676375.

PubMed Abstract | Google Scholar