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

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:

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.

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:

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.

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 (e1,e2=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:

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

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

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

which defines the forward mapping (u,v)↦(C~,Ĩ,Ṽ) with

Solving for (u, v) yields

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

To determine e₁ explicitly, we solve

with J(E₊) given in Equation 4. Setting e1=(l1,l2,1)⊤ gives the linear system

From the third row we obtain

and substituting into the first row gives

Thus,

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

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

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.

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

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

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

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:

with

for some coefficients p_i, q_i ∈ ℝ.

Expanding u², v², and uv in powers of z = u + iv and z̄=u-iv gives

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

It remains to determine the coefficients p_i and q_i at the Hopf point. By construction, u=12Ṽ, and from Equation 2 the V˙-equation is purely linear,

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

Substituting c₃ = γ/θ and c₄ = ω_0H/θ yields u˙=-ω0Hv, confirming that

To find the coefficients q₁, q₂, q₃, we consider the parametrization of e1=(c1+ic2,c3+ic4,1)⊤. Recall from Equation 11 that

and from Equation 12,

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

Because c₁, c₃, δ are constants,

Expressing C~2, C~Ĩ, and C~Ṽ in terms of (u, v) and simplifying yields

Substituting Equation 30 into Equation 29 and matching terms with v˙quad=q1u2+q2uv+q3v2 yields

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.

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

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

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

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

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

Substituting Equations 36–38 into Equation 35 yields

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

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

Equation 40 becomes

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

Thus, by choosing

we cancel all quadratic terms and obtain the cubic normal form (Equation 34). At θ = θ_H we have μ(θ_H) = iω_0H ≠ 0 and 2μ̄(θH)-μ(θH)=-3iω0H≠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:

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,

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

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

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

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

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

By setting H₃₀ = G₃₀/2μ, H12=G12/2μ̄, and H03=G03/(3μ̄-μ), we eliminate all cubic terms except w2w̄. The coefficient of the w2w̄ term is 12(G21-(μ+μ̄)H21). This term could be removed if H21=G21/(μ+μ̄), but since μ(θH)+μ̄(θH)=0 at the bifurcation point, the divisor vanishes. For this reason, w2w̄ 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:

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]:

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

which gives:

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].

We now reinterpret the Stuart–Landau normal form

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

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

For a supercritical Hopf bifurcation (a₁ > 0), the amplitude equation admits a stable limit cycle of radius r*2=λ0/a1, oscillating at frequency Ω(r*)=ω0-b1r*2.

As θ crosses θ_H, λ₀ changes sign and r^* grows from 0 for θ=θ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):

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

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

Using α2=λ0/a1 and substituting into Equation 55 yields the normal form

where

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.

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

Its oscillation frequency,

then gives

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:

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.

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.

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.

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.

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.

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.

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:

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.