Dynamics of Kawasaki Disease Pathogenesis under Stochastic Perturbations and Time-Delay Effects

1  Introduction

Kawasaki disease is an acute systemic vasculitis that predominantly occurs in young children, and may result in complications of the coronary arteries if not diagnosed and treated early. Early clinical research identified fever, mucocutaneous manifestations, lymphadenopathy and coronary artery lesions as the primary characteristics of Kawasaki disease [1]. Subsequent epidemiological studies revealed that Kawasaki disease is highly age, seasonal and geographically biased, with the greatest disease burden in East Asia [2]. A number of studies have suggested that Kawasaki disease is triggered by infectious or environmental factors in genetically predisposed individuals. While a specific pathogen has not been identified, viral, bacterial and immune mechanisms have been frequently proposed [3]. Immunological research also suggests cytokines, chemokines, endothelial dysfunction and inflammatory mediators are key factors in disease development [4]. Kawasaki disease’s inflammatory processes have also been associated with vascular endothelial damage. Increased expression of tumor necrosis factor, interleukins, vascular endothelial growth factors, adhesion molecules and chemokines have been linked to endothelial dysfunction and coronary artery disease in Kawasaki disease [5]. This suggests that a mathematical model that incorporates interactions between endothelial cells, growth factors, chemokines, and cytokines would be relevant. Mathematical models are increasingly used for investigating nonlinear mechanisms of disease. Qiang et al. constructed a differential equation model for Kawasaki disease by considering endothelial cells, vascular endothelial growth factors, adhesion factors, chemokines, and inflammatory factors. Their model exhibited rich dynamics, such as forward and backward bifurcations [6]. Guo et al. also examined the global dynamics of a Kawasaki disease model and confirmed the role of nonlinear feedbacks in sustaining the disease [7]. Delay differential equations are important for modeling non-instantaneous biological processes. In Kawasaki disease, activation of the immune system, release of cytokines, stimulation of the endothelium, and inflammatory feedback processes take time to process. The theory of delay differential equations is an appropriate mathematical formalism to describe such processes [8]. Hence, delay terms can be used to increase biological realism of Kawasaki disease models, by accounting for the time lag between inflammatory stimulation and pathological manifestation. Stochastic analysis is also important due to random variability. Immune system variability, exposure to environmental factors, disease variation and measurement uncertainty can influence disease dynamics. Mao’s theory of stochastic differential equations offers standard techniques for establishing positivity, boundedness, extinction and persistence of random dynamical systems [9]. Allen and Gray et al. demonstrated that stochastic epidemic models may lead to extinction or persistence that is not apparent in deterministic models [10,11]. Recent stochastic epidemic research has continued to employ stochastic threshold parameters to determine whether diseases go extinct or persist in the long term [12]. The incidence of Kawasaki disease was altered during the COVID-19 pandemic in various countries. A nationwide survey in Japan showed a marked decrease in Kawasaki disease incidence during the pandemic, followed by a resurgence after the pandemic ended [13–15]. Likewise, Canadian data indicated a decline in 2020–2021 and a return to baseline trends in the subsequent years [16]. This evidence indicates the possibility of the effects of environmental exposure, social contacts and infection control on Kawasaki disease incidence. The pandemic also posed diagnostic challenges because Kawasaki disease has several overlapping inflammatory symptoms with multisystem inflammatory syndrome in children. Comparative studies between pre-pandemic and pandemic periods show variations in clinical and laboratory features and treatment outcomes [17]. These observations give further impetus to use stochastic and delay models because they capture randomness, external disturbances, and time-dependent biological processes. Finite difference methods are needed to solve nonlinear stochastic delay models. Conventional finite difference methods may not always preserve positivity, boundedness and stability. Thus, nonstandard finite difference methods, initially designed to preserve qualitative features of differential equation models, are helpful in biology [18]. In stochastic models, nonstandard schemes maintain non-negativity and stability despite random disturbances, and can be applied to model stochastic Kawasaki disease dynamics. It is important to clarify the biological interpretation of the proposed model. System (1)–(4) is not a population-level epidemiological transmission model. Kawasaki disease is not modeled here as an infectious disease spreading between individuals. Instead, the model describes within-host, lesion-level interactions among endothelial cells, vascular endothelial growth factors, adhesion molecules/chemokines, and inflammatory cytokines. Therefore, the threshold quantity used in this work should be interpreted as an inflammatory feedback threshold rather than a classical epidemiological inflammatory feedback threshold.

This paper is structured as follows. Section 2 proposes the model with time delay. Section 3 provides a qualitative analysis, including the positivity and boundedness of solutions. Section 4 derives the equilibrium points and the inflammatory feedback threshold. Section 5 deals with the stability analysis of the equilibria. In Section 6, real data from 2020–2025 are analyzed to support the model. Section 7 introduces the stochastic model and its dynamical properties. Section 8 develops the stochastic NSFD scheme. Section 9 illustrates the results through graphical simulations, and Section 10 concludes the study.

2  Model Formulation

The complex interactions between endothelial cells, biochemical mediators, and inflammatory responses play an important role in Kawasaki disease pathogenesis [6]. The deterministic interaction structure used here is based on the Kawasaki disease pathogenesis model proposed by Qiang et al. [6].

In the present study, this baseline framework is extended by introducing a discrete delay in the inflammatory cytokine response, stochastic perturbations, and a structure-preserving numerical approximation. Thus, the biological interaction structure follows the previous Kawasaki disease model, while the delayed stochastic extension and the stochastic NSFD scheme are developed in this work.

It is important to emphasize that the proposed model is not a population-level epidemiological transmission model. Kawasaki disease is not modeled here as a contagious disease spreading between individuals. Instead, the model describes lesion-level inflammatory interaction mechanisms involving endothelial cells, VEGF, adhesion molecules/chemokines, and inflammatory cytokines.

The model does not explicitly divide Kawasaki disease into acute, subacute, and convalescent clinical phases. Coronary artery inflammation is also not included as a separate state variable. Instead, endothelial-cell dysfunction and inflammatory cytokine activity are used as indirect mechanistic indicators of vascular inflammatory progression.

The flow of the state variables and their interactions within the model framework are illustrated in Fig. 1.

Figure 1: Schematic representation of the proposed lesion-level inflammatory interaction model showing the relationships among endothelial cells, VEGF, adhesion molecules/chemokines, and inflammatory cytokines.

Let E(t), V(t), C(t), and P(t) denote the state variables at time t. Their biological meanings are summarized in Table 1.

The variables represent biological concentrations or activities in the lesion-level inflammatory environment of Kawasaki disease. They do not represent human population classes such as susceptible, infected, or recovered individuals. Therefore, the model describes inflammatory mechanism dynamics rather than epidemiological transmission dynamics. The interaction structure is interpreted as follows. Inflammatory cytokine activity damages healthy endothelial-cell activity through the term k1E(t)P(t). Endothelial-inflammatory interaction promotes VEGF production through the term k2E(t)P(t). The same interaction contributes to adhesion molecule/chemokine activation through k3E(t)P(t), while VEGF further stimulates adhesion molecule/chemokine activity through k4V(t). Finally, adhesion molecule/chemokine activity drives inflammatory cytokine activation after a delay through the term k5C(t−τ)e−d3τ. The variables are aggregated quantities and do not distinguish individual cytokines such as TNF-α, IL-1, IL-6, or specific adhesion molecules. The delay τ is introduced only in the inflammatory cytokine equation. It represents the time required for adhesion molecule/chemokine activity to produce a measurable inflammatory cytokine response. Biologically, this delay may include intermediate processes such as immune-cell recruitment, intracellular signaling, and cytokine amplification. Other possible delays, such as the delay from endothelial injury to VEGF production or from VEGF activity to adhesion molecule upregulation, are not modeled separately in this formulation. By introducing a discrete delay τ>0, the deterministic delayed system is given by

dE(t)dt=r+k6V(t)E(t)1+V(t)−k1E(t)P(t)−d1E(t),(1)

dV(t)dt=k2E(t)P(t)−d2V(t),(2)

dC(t)dt=k3E(t)P(t)+k4V(t)−d3C(t),(3)

dP(t)dt=k5C(t−τ)e−d3τ−d4P(t).(4)

The parameters of system (1)–(4) are summarized in

The numerical values in Table 2 are used only for qualitative simulation experiments. Where direct clinical estimates are unavailable, values are selected from biologically reasonable ranges or adapted from the baseline deterministic model of Qiang et al. [6]. Thus, the simulations should be interpreted as qualitative numerical illustrations rather than patient-specific parameter estimation.

In this formulation, r represents the baseline regeneration rate of endothelial-cell activity. The coefficients d1,d2,d3, and d4 describe the natural decay or clearance rates of the corresponding biological components. The parameter k1 represents the damaging effect of inflammatory cytokine activity on endothelial-cell activity. The parameters k2 and k3 describe the production of VEGF and adhesion molecule/chemokine activity due to endothelial-inflammatory interaction, while k4 describes the additional stimulation of adhesion molecule/chemokine activity by VEGF. The parameter k5 controls delayed inflammatory cytokine activation, and k6 describes VEGF-mediated endothelial feedback. The delay parameter τ>0 denotes the average time between adhesion molecule/chemokine activation and the subsequent inflammatory cytokine response. Thus, the delayed term C(t−τ) represents earlier adhesion molecule/chemokine activity that contributes to cytokine activation at time t. The initial conditions are prescribed on the interval [−τ,0] as

E(θ)=ϕ1(θ),V(θ)=ϕ2(θ),C(θ)=ϕ3(θ),P(θ)=ϕ4(θ),−τ≤θ≤0,(5)

where

ϕi(θ)≥0,i=1,2,3,4.

3  Qualitative Analysis

In this section, we study the positivity and boundedness of system (1)–(4). Since all state variables represent biological concentrations, the solutions must remain nonnegative and bounded for the model to be biologically meaningful.

R+4={(E,V,C,P)∈R4:E≥0, V≥0, C≥0, P≥0}.

For the delayed system, we consider the phase space

𝒞+=C([−τ,0],R+4),

with nonnegative initial functions given by (5). In the this manuscript, we do not define the feasible region by the constants ME,MV,MC, and MP, because such bounds are mutually dependent and therefore circular. Instead, we construct an absorbing positively invariant set using a Lyapunov-type functional.

Theorem 1 (Positivity of solutions): For any nonnegative initial functions

ϕi(θ)≥0,−τ≤θ≤0,i=1,2,3,4,

the solution (E(t),V(t),C(t),P(t)) of system (1)–(4) remains nonnegative for all t≥0.

Proof: We verify the vector field on the boundary of the nonnegative orthant. If E(t)=0, then from (1),

dE(t)dt|E=0=r>0.

Thus, E(t) cannot become negative. If V(t)=0, then

dV(t)dt|V=0=k2E(t)P(t)≥0.

If C(t)=0, then

dC(t)dt|C=0=k3E(t)P(t)+k4V(t)≥0.

Finally, if P(t)=0, then

dP(t)dt|P=0=k5C(t−τ)e−d3τ≥0.

Therefore, the vector field points inward or is tangent on each boundary component of R+4. Hence, every solution starting from 𝒞+ remains in R+4 for all t≥0. □

Theorem 2 (Boundedness and positive invariance): Assume that d1>k6. Then every nonnegative solution of system (1)–(4) is ultimately bounded. Moreover, the system admits an absorbing positively invariant set in 𝒞+.

Proof: Since

V(t)1+V(t)≤1,

the first equation gives

dE(t)dt≤r−(d1−k6)E(t).

Because d1>k6, comparison theory gives

lim supt→∞E(t)≤rd1−k6.

Thus, E(t) is ultimately bounded.

To obtain a non-circular bound for all variables, define

𝒲(t)=E(t)+aV(t)+bC(t)+cP(t)+q∫t−τte−μ(t−s)C(s)ds,

where a,b,c,q,μ>0. Differentiating 𝒲(t) along solutions gives

𝒲˙(t)=r+k6V(t)E(t)1+V(t)−k1E(t)P(t)−d1E(t)+a(k2E(t)P(t)−d2V(t))+b(k3E(t)P(t)+k4V(t)−d3C(t))+c(k5e−d3τC(t−τ)−d4P(t))+qC(t)−qe−μτC(t−τ)−μq∫t−τte−μ(t−s)C(s)ds.

Choose

q=ck5e−d3τeμτ.

Then the delayed term C(t−τ) is cancelled. Using

k6V(t)E(t)1+V(t)≤k6E(t),

we obtain

𝒲˙(t)≤r−(d1−k6)E(t)−(ad2−bk4)V(t)−(bd3−q)C(t)−cd4P(t)−(k1−ak2−bk3)E(t)P(t)−μq∫t−τte−μ(t−s)C(s)ds.

Choose a,b,c,μ>0 such that

ak2+bk3<k1,bk4<ad2,ck5e−d3τeμτ<bd3.

Then all negative coefficients above are well-defined and positive. Hence, there exists η>0 such that

𝒲˙(t)≤r−η𝒲(t).

By comparison,

𝒲(t)≤𝒲(0)e−ηt+rη(1−e−ηt),

and therefore

lim supt→∞𝒲(t)≤rη.

Since 𝒲(t) is positive and controls E(t),V(t),C(t), and P(t), all state variables are ultimately bounded.

For any R≥r/η, define

ΩR={φ∈𝒞+:𝒲(φ)≤R},

where

𝒲(φ)=φ1(0)+aφ2(0)+bφ3(0)+cφ4(0)+q∫−τ0eμθφ3(θ)dθ.

On the boundary 𝒲=R,

𝒲˙(t)≤r−ηR≤0.

Therefore, the vector field points inward on the boundary of ΩR. Thus, ΩR is positively invariant. Since every solution eventually enters such a set, ΩR is also absorbing. □

4  Equilibrium Points and Inflammatory Feedback Threshold

Since the delay is constant and equilibrium states are time independent, we have C(t−τ)=C∗ at equilibrium. Therefore, the introduction of the discrete delay does not change the algebraic procedure used to determine the equilibrium points. The equilibrium structure follows the corresponding analysis of Qiang et al. [6]. In the present work, we recall the equilibrium expressions only for completeness and to establish notation for the subsequent delay-dependent stability analysis, stochastic formulation, and numerical approximation. The main modification introduced by the delay appears through the factor e−d3τ, which changes the effective threshold quantity. At equilibrium, all time-dependent variables become constant. Hence, the delay term satisfies C(t−τ)=C∗. From (1)–(4), the inflammatory-factor-free equilibrium is

ℰ0=(E0,V0,C0,P0)=(rd1,0,0,0).(6)

Following the threshold construction used for the corresponding Kawasaki disease interaction model, we define the inflammatory feedback threshold quantity by considering the coupled variables V, C, and P. This quantity is not a population-level epidemiological reproduction number. Instead, it measures the ability of inflammatory feedback among VEGF, adhesion molecules/chemokines, and inflammatory cytokines to sustain a positive inflammatory state. The effective inflammatory feedback threshold is obtained by applying the next-generation matrix method to the variables V,C, and P. At ℰ0, we have

F=(00k2E000k3E0000),𝒱=(d200−k4d300−k5e−d3τd4).

Therefore, the inflammatory feedback threshold is

ℛ0=ρ(F𝒱−1)=rk5e−d3τ(k2k4+d2k3)d1d2d3d4.(7)

Thus, a larger value of this threshold indicates stronger inflammatory feedback. The factor e−d3τ shows that the delayed chemokine-mediated response modifies the effective inflammatory feedback strength. The factor e−d3τ shows that the delay reduces the effective contribution of chemokine-mediated inflammatory activation. Thus, a larger delay τ decreases ℛ0.

For a positive inflammatory equilibrium

ℰ∗=(E∗,V∗,C∗,P∗),

where V∗,C∗,P∗>0, the equilibrium relations give

E∗=d2d3d4ed3τk5(k2k4+d2k3)=rd1ℛ0.(8)

Moreover,

P∗=d2V∗k2E∗,C∗=d4ed3τk5P∗.(9)

The value of V∗ is obtained from the quadratic equation

A(V∗)2+BV∗+D=0,(10)

where

A=k1d2k2,B=k1d2k2+d1E∗−r−k6E∗,D=d1E∗−r.

Hence,

V∗=−B+B2−4AD2A,(11)

provided that B2−4AD≥0 and V∗>0. Once V∗ is determined, P∗ and C∗ follow directly from (9). By using (8),

E∗=rd1ℛ0,

the coefficients in (10) become

A=k1d2k2,

B=k1d2k2+rℛ0−r−k6rd1ℛ0,

and

D=rℛ0−r=r(1−ℛ0)ℛ0.

Therefore,

V∗=−(k1d2k2+rℛ0−r−k6rd1ℛ0)+(k1d2k2+rℛ0−r−k6rd1ℛ0)2−4(k1d2k2)(r(1−ℛ0)ℛ0)2(k1d2k2).(12)

Equivalently,

V∗=k22k1d2[−B+B2−4k1d2r(1−ℛ0)k2ℛ0],(13)

where

B=k1d2k2+rℛ0−r−k6rd1ℛ0.

After obtaining V∗, the remaining components of the inflammatory equilibrium are

P∗=d2V∗k2E∗=d1d2ℛ0rk2V∗,(14)

and

C∗=d4ed3τk5P∗=d1d2d4ℛ0ed3τrk2k5V∗.(15)

Hence, the inflammatory equilibrium is

ℰ∗=(rd1ℛ0, V∗, d1d2d4ℛ0ed3τrk2k5V∗, d1d2ℛ0rk2V∗).(16)

5  Stability Analysis

In this section, we analyze the local and global stability of the inflammatory-factor-free and inflammatory equilibria of the system (1)–(4).

Theorem 3: The inflammatory-factor-free equilibrium

ℰ0=(rd1,0,0,0)

of system (1)–(4) is locally asymptotically stable if ℛ0<1, and unstable if ℛ0>1.

Proof: The local stability of ℰ0 is determined from the linearized delay system. At

ℰ0=(rd1,0,0,0),

the characteristic equation is

(λ+d1)[(λ+d2)(λ+d3)(λ+d4)−k5E0e−d3τ(k2k4+k3(λ+d2))e−λτ]=0.

Thus, one characteristic root is λ=−d1<0. The remaining roots are given by

(λ+d2)(λ+d3)(λ+d4)=k5E0e−d3τ(k2k4+k3(λ+d2))e−λτ.

Let

K=k5E0e−d3τ.

Then the above equation becomes

(λ+d2)(λ+d3)(λ+d4)=K(k2k4+k3(λ+d2))e−λτ.

We first prove stability for ℛ0<1. Suppose, to the contrary, that there exists a characteristic root λ with Re⁡(λ)≥0. Taking moduli on both sides gives

|λ+d2||λ+d3||λ+d4|=K|k2k4+k3(λ+d2)||e−λτ|.

Since Re⁡(λ)≥0, we have

|e−λτ|≤1,|λ+di|≥di,i=2,3,4.

Therefore,

|λ+d2||λ+d3||λ+d4|≥|λ+d2|d3d4.

On the other hand,

K|k2k4+k3(λ+d2)||e−λτ|≤K(k2k4+k3|λ+d2|).

Hence, a root with Re⁡(λ)≥0 would require

|λ+d2|d3d4≤K(k2k4+k3|λ+d2|).

Let x=|λ+d2|. Since x≥d2, this inequality becomes

xd3d4≤K(k2k4+k3x).

Equivalently,

x(d3d4−Kk3)≤Kk2k4.

However, from ℛ0<1, we have

K(k2k4+d2k3)d2d3d4<1,

which implies

K(k2k4+d2k3)<d2d3d4.

Hence,

d2(d3d4−Kk3)>Kk2k4.

Since x≥d2, it follows that

x(d3d4−Kk3)≥d2(d3d4−Kk3)>Kk2k4.

This contradicts the necessary inequality

x(d3d4−Kk3)≤Kk2k4.

Therefore, no characteristic root can satisfy Re⁡(λ)≥0. Hence, all characteristic roots have negative real parts, and ℰ0 is locally asymptotically stable when ℛ0<1. If ℛ0>1, then

K(k2k4+d2k3)>d2d3d4.

Define

G(λ)=(λ+d2)(λ+d3)(λ+d4)−K(k2k4+k3(λ+d2))e−λτ.

Then

G(0)=d2d3d4−K(k2k4+d2k3)<0.

Moreover,

limλ→+∞G(λ)=+∞.

By continuity, there exists a positive real root λ>0. Therefore, the inflammatory-factor-free equilibrium ℰ0 is unstable when ℛ0>1. □

Theorem 4: Let ℰ∗ be a positive inflammatory equilibrium of system (1)–(4). If all characteristic roots of the linearized delay system at ℰ∗ have strictly negative real parts, then ℰ∗ is locally asymptotically stable. If characteristic roots lie on the imaginary axis, the linearized analysis is inconclusive and further center-manifold or bifurcation analysis is required.

Proof: The local behavior of system (1)–(4) near the positive inflammatory equilibrium ℰ∗=(E∗,V∗,C∗,P∗) is determined by the Jacobian matrix evaluated at ℰ∗. Since the last equation contains the delayed term C(t−τ), the linearization introduces the factor e−λτ. Therefore, the delay-Jacobian matrix at ℰ∗ is

J(ℰ∗)=(a11a120a14a21−d20a24a31k4−d3a3400k5e−d3τe−λτ−d4),

where

a11=k6V∗1+V∗−k1P∗−d1,a12=k6E∗(1+V∗)2,a14=−k1E∗,

a21=k2P∗,a24=k2E∗,a31=k3P∗,a34=k3E∗.

The eigenvalues are obtained from the characteristic equation

det(λI−J(ℰ∗))=0.

Thus,

|λ−a11−a120−a14−a21λ+d20−a24−a31−k4λ+d3−a3400−k5e−d3τe−λτλ+d4|=0.(17)

Expanding (17) along the fourth row gives

(λ+d4)|λ−a11−a120−a21λ+d20−a31−k4λ+d3|−k5e−d3τe−λτ|λ−a11−a12−a14−a21λ+d2−a24−a31−k4−a34|=0.(18)

The first determinant in (18) is

|λ−a11−a120−a21λ+d20−a31−k4λ+d3|=(λ+d3)[(λ−a11)(λ+d2)−a12a21].

The second determinant is

|λ−a11−a12−a14−a21λ+d2−a24−a31−k4−a34|=a14[a21k4+a31(λ+d2)]+a24[k4(λ−a11)+a12a31]+a34[(λ−a11)(λ+d2)−a12a21].

Hence, the characteristic equation can be written as

(λ+d4)(λ+d3)[(λ−a11)(λ+d2)−a12a21]−k5e−d3τe−λτQ(λ)=0,(19)

where

Q(λ)= a14[a21k4+a31(λ+d2)]+a24[k4(λ−a11)+a12a31]+a34[(λ−a11)(λ+d2)−a12a21].(20)

Now, expanding the polynomial part, we obtain

(λ−a11)(λ+d2)−a12a21=λ2+(d2−a11)λ−(a11d2+a12a21).

Therefore,

(λ+d4)(λ+d3)=λ2+(d3+d4)λ+d3d4.

Thus, the non-delay polynomial part becomes

P(λ)=λ4+A1λ3+A2λ2+A3λ+A4,(21)

where

A1=d2+d3+d4−a11,

A2=d3d4+(d3+d4)(d2−a11)−(a11d2+a12a21),

A3=d3d4(d2−a11)−(d3+d4)(a11d2+a12a21),

and

A4=−d3d4(a11d2+a12a21).

We emphasize that the above spectral condition is used only as a sufficient mathematical condition for local asymptotic stability of the linearized delayed system. It is not used as a direct biological criterion. Biologically, local stability means that small perturbations in endothelial-cell activity, VEGF, adhesion molecules/chemokines, and inflammatory cytokines decay back toward the steady inflammatory state. Similarly, the delay-dependent polynomial Q(λ) can be written as

Q(λ)=B2λ2+B1λ+B0,(22)

where

B2=a34,

B1=a14a31+a24k4+a34(d2−a11),

and

B0=a14(a21k4+a31d2)+a24(−k4a11+a12a31)−a34(a11d2+a12a21).

Consequently, the characteristic equation takes the polynomial-delay form

λ4+A1λ3+A2λ2+A3λ+A4−k5e−d3τe−λτ(B2λ2+B1λ+B0)=0.(23)

The local stability of ℰ∗ is determined by the roots of (23). If all roots λ satisfy

Re⁡(λ)<0,

then every sufficiently small perturbation around ℰ∗ decays to zero as t→∞. Hence, the positive inflammatory equilibrium is locally asymptotically stable.

Furthermore, for λ=iω, we have

|e−iωτ|=1.

Hence,

|k5e−d3τe−iωτ(B2(iω)2+B1(iω)+B0)|=k5e−d3τ|B2(iω)2+B1(iω)+B0|.

Therefore, a sufficient stability condition is

k5e−d3τ|B2(iω)2+B1(iω)+B0|<|(iω)4+A1(iω)3+A2(iω)2+A3(iω)+A4|,∀ω≥0.

Under this condition, the delayed feedback term is dominated by the stable polynomial part. Thus, no characteristic root enters the right half-plane, and ℰ∗ remains locally asymptotically stable. Conversely, if at least one characteristic root has positive real part, then some perturbations grow near ℰ∗, and the positive inflammatory equilibrium is unstable. □

Theorem 5: The inflammatory-factor-free equilibrium

ℰ0=(rd1,0,0,0)

of system (1)–(4) is locally asymptotically stable if ℛ0<1, and unstable if ℛ0>1, provided that no characteristic root lies on the imaginary axis.

Proof: The local stability of ℰ0 is determined from the linearized delay system. At

ℰ0=(rd1,0,0,0),

the characteristic equation is

(λ+d1)[(λ+d2)(λ+d3)(λ+d4)−k5ℰ0e−d3τ(k2k4+k3(λ+d2))e−λτ]=0.

Hence, one characteristic root is λ=−d1<0. The remaining characteristic roots satisfy

(λ+d2)(λ+d3)(λ+d4)−k5ℰ00e−d3τ(k2k4+k3(λ+d2))e−λτ=0.

At λ=0, the left-hand side becomes

d2d3d4−k5ℰ0e−d3τ(k2k4+d2k3).

Using ℰ0=r/d1, we obtain

d2d3d4−k5ℰ0e−d3τ(k2k4+d2k3)=d2d3d4(1−ℛ0).

Therefore, ℛ0 acts as a local threshold quantity near ℰ0. If ℛ0>1, the characteristic equation admits a positive real root, and hence ℰ0 is unstable. If ℛ0<1, and all characteristic roots of the linearized delayed system have negative real parts, then ℰ0 is locally asymptotically stable by the local stability theory of delay differential equations. Because the underlying Kawasaki disease interaction model may exhibit backward bifurcation, the condition ℛ0<1 is not sufficient by itself to guarantee global asymptotic stability. Thus, in the revised manuscript, we do not claim global convergence to ℰ0 based only on ℛ0<1. □

Remark 1: The above result is a local stability result. In the present model, ℛ0<1 should not be interpreted as a sufficient condition for global elimination of inflammatory activity. Since backward bifurcation may occur, a positive inflammatory equilibrium can coexist with the inflammatory-factor-free equilibrium for some parameter values. Therefore, additional restrictions are required before any global stability conclusion can be made.

Theorem 6: Let

ℰ∗=(E∗,V∗,C∗,P∗)

be a positive inflammatory equilibrium of system (1)–(4). If all characteristic roots of the linearized delay system at ℰ∗ have strictly negative real parts, then ℰ∗ is locally asymptotically stable.

Proof: Let

u1(t)=E(t)−E∗,u2(t)=V(t)−V∗,u3(t)=C(t)−C∗,u4(t)=P(t)−P∗

be small perturbations around the positive inflammatory equilibrium ℰ∗. The local behavior of system (1)–(4) near ℰ∗ is determined by the linearized delayed system. Since the fourth equation contains the delayed term C(t−τ), the linearization introduces the exponential factor e−λτ in the characteristic equation. The corresponding characteristic equation can be written as

λ4+A1λ3+A2λ2+A3λ+A4−k5e−d3τe−λτ(B2λ2+B1λ+B0)=0,

where the coefficients Ai and Bi are obtained from the Jacobian matrix evaluated at ℰ∗.

If every characteristic root λ of this equation satisfies

Re⁡(λ)<0,

then the zero solution of the linearized perturbation system is asymptotically stable. Therefore, by the local stability theory of delay differential equations, the positive inflammatory equilibrium ℰ∗ is locally asymptotically stable. We emphasize that this result is local. It does not imply global asymptotic stability in the whole feasible region Ω. A global Lyapunov result would require additional parameter restrictions ensuring that all nonlinear cross-interaction terms are controlled. Since such restrictions are not derived here, no global stability claim is made for ℰ∗. □

6  Reported Incidence Data as Qualitative Motivation

In this section, we summarize recently reported incidence data for Kawasaki disease during and after the COVID-19 pandemic. These data are used only as qualitative motivation for considering variability and delayed biological responses in Kawasaki disease. They are not used to calibrate, estimate, or statistically validate the proposed stochastic delayed model. In particular, the incidence data are population-level observations, whereas system (1)–(4) describes lesion-level inflammatory interactions among endothelial cells, VEGF, adhesion molecules/chemokines, and inflammatory cytokines.

The summary focuses on reported incidence rates among children aged 0 to 4 years, because this age group has the highest reported burden of Kawasaki disease. Published studies indicate that Kawasaki disease incidence decreased during the early COVID-19 pandemic period and subsequently increased after the relaxation of public health measures. In Canada, the pre-pandemic incidence was approximately 22.5 cases per 100,000 children aged 0–4 years, while it declined to 15.8 cases per 100,000 during 2020–2021. The incidence later returned toward the pre-pandemic pattern, reaching approximately 24.1 cases per 100,000 during 2023–2024 [16] (See Table 3). Similarly, data from the national Japanese survey showed a decrease during the pandemic, followed by a marked resurgence after the relaxation of public health measures [13,15].

The percentage change in incidence was calculated using

%Δ=I2−I1I1×100,

where I1 and I2 denote the baseline and comparison incidence rates, respectively. For Canada, the incidence decreased from 22.5 to 15.8 per 100,000, corresponding to

15.8−22.522.5×100=−29.78%.

This represents an approximately 30% reduction during the early pandemic period. In contrast, the incidence increased from 15.8 in 2020–2021 to 24.1 in 2023–2024, giving

24.1−15.815.8×100=52.53%.

Thus, the reported post-pandemic period showed a rebound in Kawasaki disease incidence.

For Japan, the reported incidence increased from 269.3 in 2021 to 426.7 per 100,000 children aged 0–4 years in 2023. The relative increase is

426.7−269.3269.3×100=58.45%.

This increase is consistent with reports of a resurgence after the relaxation of pandemic-related restrictions. However, these population-level incidence patterns should not be interpreted as direct validation of the lesion-level mathematical model.

Figs. 2 and 3 show regional and temporal differences in reported Kawasaki disease incidence. The decline during 2020–2021 may be associated with reduced exposure to infectious or environmental triggers during the COVID-19 pandemic, while the later increase may reflect changes following the relaxation of public health restrictions. Japan shows a stronger reported increase than Canada, suggesting region-specific incidence patterns.

Figure 2: Reported temporal trends of Kawasaki disease incidence in Canada and Japan. The figure provides qualitative motivation only and is not used for model validation.

Figure 3: Heatmap representation of reported Kawasaki disease incidence intensity across regions and years. The figure illustrates regional and temporal variation only.

These figures should be interpreted only as qualitative background motivation. They do not verify the theoretical stability results, do not support a stochastic threshold, and do not validate the stochastic NSFD scheme. The purpose of including these data is to show that Kawasaki disease incidence can vary substantially across time and region, which motivates the mathematical consideration of variability and delayed biological response mechanisms.

7  Stochastic Model

The stochastic perturbations are introduced to represent random variability in the biological activity of each model component. They are not intended to model a specific infectious agent or environmental exposure as a separate dynamical variable. Instead, the multiplicative noise terms describe unresolved fluctuations in endothelial response, VEGF activity, adhesion molecule/chemokine activation, and inflammatory cytokine activity. Such fluctuations may reflect patient-specific immune variability, local microenvironmental differences, measurement uncertainty, and unmodeled triggering events such as infectious or environmental exposures. To incorporate random fluctuations arising from biological variability, environmental uncertainty, and measurement noise, we extend the deterministic delay system (1)–(4) into a stochastic delay differential equation model. Let σ1,σ2,σ3,σ4≥0 denote the noise intensities associated with E(t),V(t),C(t), and P(t), respectively. Let Bi(t), i=1,2,3,4, be mutually independent standard Brownian motions. The stochastic delay model is given by

dE(t)=[r+k6V(t)E(t)1+V(t)−k1E(t)P(t)−d1E(t)]dt+σ1E(t)dB1(t),(24)

dV(t)=[k2E(t)P(t)−d2V(t)]dt+σ2V(t)dB2(t),(25)

dC(t)=[k3E(t)P(t)+k4V(t)−d3C(t)]dt+σ3C(t)dB3(t),(26)

dP(t)=[k5C(t−τ)e−d3τ−d4P(t)]dt+σ4P(t)dB4(t).(27)

Here, the stochastic perturbation terms are chosen in proportional form. Thus, the magnitude of the random fluctuation depends on the current level of the corresponding biological variable. In particular, σ1E(t)dB1(t) represents random variability in endothelial-cell activity and local injury/repair responses; σ2V(t)dB2(t) represents fluctuations in VEGF activity; σ3C(t)dB3(t) represents variability in adhesion molecule/chemokine activation; and σ4P(t)dB4(t) represents random fluctuations in inflammatory cytokine activity. The parameters σi quantify the intensity of these fluctuations.

7.1 Positivity and Boundedness of the Stochastic Delayed Model

Let (Ω,ℱ,P) be a complete probability space equipped with a filtration {ℱt}t≥0 satisfying the usual conditions. Consider the stochastic delayed model (24)–(27) defined on this space.

Let

X(t)=(E(t),V(t),C(t),P(t)),

and define the norm

|X(t)|=E2(t)+V2(t)+C2(t)+P2(t).

We consider the class C2,1(R+4×(0,∞);R+) of all nonnegative functions that are twice continuously differentiable in X and once in t.

The stochastic delay system can be written in compact form as

dX(t)=F(X(t),t)dt+G(X(t),t)dW(t),(28)

where F(X,t) and G(X,t) denote the drift and diffusion terms, respectively, and W(t) is a standard Brownian motion vector.

Define the differential operator ℒ associated with (28) by

ℒV(X,t)=Vt+∑i=14Fi(X,t)∂V∂Xi+12∑i,j=14(GGT)ij∂2V∂Xi∂Xj.

Theorem 7 (Positivity and global existence of the stochastic delayed model): Assume that the initial functions satisfy

E(θ)>0,V(θ)>0,C(θ)>0,P(θ)>0,−τ≤θ≤0.

Suppose that d1>k6 and that there exist positive constants a,b,c,μ such that

ak2+bk3<k1,bk4<ad2,ck5e−d3τeμτ<bd3.

Then system (24)–(27) admits a unique global positive solution. Moreover,

(E(t),V(t),C(t),P(t))∈R+4for all t≥0a.s.

Proof: The drift and diffusion coefficients of (24)–(27) are locally Lipschitz on R+4. Hence, for any positive initial function, there exists a unique local solution on [0,τe), where τe denotes the explosion time. We emphasize that the global linear growth condition is not assumed here, because the drift contains bilinear terms such as E(t)P(t). We first prove positivity. Up to the explosion time, the first equation can be written in the linear form

dE(t)=[r+A1(t)E(t)]dt+σ1E(t)dB1(t),

where

A1(t)=k6V(t)1+V(t)−k1P(t)−d1.

By the variation-of-constants formula,

E(t)=Φ1(t)[E(0)+∫0tΦ1−1(s)rds],

where

Φ1(t)=exp⁡(∫0t(A1(s)−σ122)ds+σ1B1(t))>0.

Therefore, E(t)>0 for t<τe. Similarly, the equation for V(t) is

dV(t)=[k2E(t)P(t)−d2V(t)]dt+σ2V(t)dB2(t).

Thus,

V(t)=Φ2(t)[V(0)+∫0tΦ2−1(s)k2E(s)P(s)ds]>0,