Introduction

The study of epidemic spread dynamics, apart from being of theoretical interest, is an important precursor for developing effective interventions for mitigation based on control-theoretic analysis. The complex dynamics of epidemics was exemplified by COVID-19, marked as it was by dynamical characteristics such as infection spikes induced by superspreader events, varying transmissibility of mutants of the causative virus SARSCOV-2, and so on. More specifically, understanding dynamic instabilities triggered by the nonlinear and stochastic characteristics of infection spread and the consequent emergence of self-organized spread patterns is the facet of the aforementioned complex dynamics that fundamentally motivates this paper. An effective approach to the study of instabilities is to consider an epidemic as a nonlinear, stochastic dynamical system within the framework of a compartmental epidemic model [1]. Here, the nonlinearity emerges in mathematical models of a susceptible population becoming infected. On the other hand, stochasticity is essential to describe uncertainties (for example, in human behavior and pathogen transmission characteristics) that critically influence the spread dynamics. Notably, describing spatiotemporal dynamics requires a coupled partial differential equation (PDE) framework for the compartmental model since ordinary differential equation (ODE) models only account for temporal dynamics. The effectiveness of the PDE approach for predicting epidemic spread was highlighted using reaction-diffusion type PDEs in our previous work [2,3]; an optimization scheme enabled the PDE framework to first 1 Corresponding Author. learn the critical model parameters from actual COVID-19 spread data for a 30 day period and then predict the spread dynamics for the succeeding 15 days. Validating the approach using COVID-19 data from the entire state of Ohio, United States, as well as the county of Hamilton, Ohio, we also found that the numerical predictions of COVID-19 spread obtained from the PDE model were more accurate for the smaller geographical region of the county than the entire state. These results also motivate the study of instabilities in PDE epidemic models reported in the present work. Instabilities in reaction-diffusion PDEs and their relationship to steady-state pattern formation were first studied by Alan Turing. He discovered that a homogeneous steady-state solution of a reactiondiffusion PDE can develop instabilities owing to the disparity of diffusion coefficients of the two reacting species considered [4,5]. Standing apart for its universality, the Turing instability has been used to explain a variety of phenomena -particularly in the biological sciences -such as the formation of striped patterns on zebras [6][7][8][9]. Remarkably, experimental evidence for the instability was reported in the 1990s [10] motivating further advances [4,11]. However, Turing's stability analysis is conspicuously based on linearized dynamics near the steady state. Therefore, despite its elegance in the case of linear reaction-diffusion PDEs, the Turing analysis is not directly applicable to nonlinear systems. Nonetheless, it has motivated higher-order perturbative analyses aimed at understanding the role of nonlinearity and noise in triggering instabilities, leading to self-organized pattern formation in nonlinear stochastic reaction-diffusion PDEs [12]. In recent work, we studied the onset of dynamic instabilities and consequent pattern formation in a stochastic PDE model of epidemic spread, invoking higher-order perturbative analysis [13,14].

The key findings were: (1) Higher-order perturbative analysis uncovered dynamic instabilities induced by stochastic environmental effects -represented by white noise -in the presence of nonlinear infection forces. Remarkably, a lower-order analysis could not reveal these instabilities. (2) Saturation parameters in nonlinear infection forces in tandem with diffusion coefficients and white noise intensity can influence the onset of instability. (3) Self-organized, steady-state patterns of infection spread were found, corresponding to the instabilities. In particular, we reported the emergence of diffusion-driven patterns in the deterministic case and white noise-induced patterns in cases when diffusion alone does not induce steady-state patterns (i.e. noise-induced patterns even in the absence of the classical Turing patterns). In this paper, we investigate the case of stochastic environmental effects on the infection spread dynamics represented by the Ornstein-Uhlenbeck (O-U) stochastic process. The fundamental point that motivates this effort is that the non-zero correlation time of Ornstein-Uhlenbeck noise renders it a significantly better representation of real-world uncertainties driving an epidemic spread than the uncorrelated white noise process. Indeed, it is wellrecognized that white noise, while mathematically tractable and attractive, is often an inadequate approximation of stochasticity in the real world [15]. It is also of interest to note that the O-U noise process is a non-Markovian process (a process with "memory"). Note that while most real-world random phenomena are more accurately modeled as non-Markovian processes [16], the Markovian approximation (with its idealized, vanishing correlation time) is often used since it more readily yields analytic results. In mathematical epidemiology, while correlated noise has been utilized in a few studies on compartmental models of epidemic spread [17,18], it is yet to be understood in the context of reaction-diffusion-type epidemic models. Here we seek to understand the fundamental role of the finite correlation time of the O-U noise process on the self-organized behavior and to analyze its effects on the emergence of stationary patterns. Furthermore, we show that our previous results corresponding to the white noise case are recovered in the limit as the correlation time of the O-U noise process tends to zero [13,14]. To reiterate, the study of instabilities in the PDE epidemic model under Ornstein-Uhlenbeck noise is the novel contribution of this work. Next, we briefly discuss the rationale for considering the nonlinearity and noise in epidemic dynamics. In compartmental models (of which our PDE model is an example), the population is partitioned into disjoint subsets of Susceptible (S), Infected (I), and Recovered (R) individuals; epidemic dynamics is characterized by the I compartment gaining individuals owing to interaction with the S compartment. This transition is mathematically represented by the infection force. Based on empirical observations of saturation effects in epidemic spread patterns, the infection force is often well represented by strongly nonlinear functions [19][20][21]. Therefore one obtains a nonlinear PDE model. On the other hand, random environmental effects are well-recognized factors driving epidemic spread [22,23]. Therefore nonlinearity and noise ought to be taken into account in stability analysis. The rest of the paper is set as follows. The PDE model and the mathematical framework for the analysis are presented in Section 2. The results are presented in Section 3. The paper concludes in Section 4 with a brief discussion of the results and an outlook for further research.

PDE Model and Mathematical Framework:

The full system of reaction-diffusion type PDEs in the fourcompartment model of epidemic dynamics studied in our previous work consists of four coupled equations. Each PDE traces the spatiotemporal evolution of the Susceptible (𝑆), Latent (𝐿), Infected (𝐼), and Recovered (𝑅) population densities corresponding to the four distinct compartments [13,14,24,25]. Here, we consider a 𝑆𝐼 𝑅 model and analyze the reduced, coupled system of the 𝑆 and 𝐼 PDEs. We note here that, in the absence of reinfections that lie out-side the scope of the present analysis (we will consider reinfection dynamics as a separate case in future work), the 𝑅 compartment only serves as a reservoir for the recovered individuals and does not significantly contribute to the overall dynamics. Therefore, studying the reduced system can yet reveal interesting aspects of the coupled 𝑆-𝐼 dynamics. We now turn to the rate of infection -the source of nonlinearity in the problem. In the simplest models, the rate of infection is defined as 𝛽(𝐼) = 𝛽 0 𝐼, where the constant 𝛽 0 represents the transmission rate, and 𝐼 is the infected population density. Consequently, the infection force -the key coupling term that describes the transition from the 𝑆 to the 𝐼 compartment -is taken to be 𝛽 0 𝑆𝐼. Note that the infection force thus defined is multiplicative; hence a nonlinear function of 𝑆 and 𝐼. However, practical considerations suggest that infection forces with growth constraints on 𝐼 owing to saturation effects can provide more realistic models of the 𝑆 to 𝐼 transition than the unbounded growth of 𝐼 produced by 𝛽 0 𝑆𝐼. These considerations include the observation that, in the absence of a correspondingly large increase in the number of susceptible individuals, a rapidly spreading infection will run out of susceptibles to infect, and saturate. Moreover, factors such as interventions by public health administrations, and human behavioral changes induced by increased awareness as the infection spreads across communities, can also inhibit and saturate the spread [26]. Accounting for such saturation effects, the infection force can be modeled as: (see, for instance, [20])

where 𝛼 is a positive constant that represents the saturation parameter.

2.1 Model:. We consider a nonlinear infection force with saturation effects (see, Eq. ( 1)) to obtain a reduced system of coupled, reaction-diffusion PDEs for 𝑆 and 𝐼 as:

𝜕𝐼 𝜕𝑡

where the birth rate is denoted by 𝑏, the death rate is 𝑑 in the susceptible and infected population (owing to reasons other than the considered epidemic), and 𝛾 is the recovery rate in the infected population (represents individuals who recover due to vaccination from the epidemic). The diffusion coefficient corresponding to the susceptible population density 𝑆 is taken to be 𝐷. Note here that, the diffusion coefficient in the PDE for 𝐼 is assumed to be unity for ease of analysis. In other words, while we thus scale the diffusion parameters, the analysis that follows will be valid for other constant diffusion coefficients, say, 𝐷 1 and 𝐷 2 for the 𝑆 and 𝐼 PDEs, respectively. The Laplacian operator is ∇ 2 , and 𝜉 (𝑡, 𝑥, 𝑦) is the spatiotemporal Ornstein-Uhlenbech (O-U) process, which at a point (𝑥, 𝑦) in space can be defined as [15] 𝜏 𝜕 𝜕𝑡 𝜉 (𝑥, 𝑦, 𝑡) = -𝜉 (𝑥, 𝑦, 𝑡) + 𝜂(𝑥, 𝑦, 𝑡),

where the autocorrelation function of Gaussian white noise 𝜂(𝑡, 𝑥, 𝑦) reads as

Here 𝜎 is the noise intensity and the 𝛿(.) represent the Dirac delta function of its argument. The correlation function of the spatiotemporal O-U process is given by: 2 / Transactions of the ASME

We note that this O-U noise term in Eq. ( 2) represents the random effects of the environment on the dynamics of 𝐼. Consider (𝑆 0 , 𝐼 0 ), the homogeneous steady state solution of the stochastic PDE system described by Eq. ( 2). The equilibrium quantities 𝑆 0 , and 𝐼 0 may be obtained by setting the right-hand side of Eq. ( 2) to zero, in the absence of diffusion and noise terms. Thus (𝑆 0 , 𝐼 0 ) simultaneously satisfies:

The nontrivial solution of Eq. ( 6) is:

It is important to note here that, in the absence of both the external noise (perturbing the PDE for 𝐼) and the diffusion term, the solution (𝑆 0 , 𝐼 0 ) is stable. However, it is to be expected that competing populations in reaction-diffusion systems, if allowed to diffuse, will self-organize leading to unique patterns. In other words, pattern formation occurs if diffusion causes instability in the homogeneous steady solution; indeed this remarkable phenomenon is termed the Turing instability. Importantly though, while diffusion alone might be insufficient for creating instabilities in a reaction-diffusion epidemic model, diffusing population densities may undergo instabilities in the presence of noise of even small intensity. This interplay between noise and diffusion becomes even more intriguing in the presence of saturation effects and here we study this with a focus on O-U noise.

Stability and Moments.

For stability analysis, we consider small perturbations to the uniform steady solutions, i.e., 𝑆 → 𝑆 0 + 𝛿𝑆, and 𝐼 → 𝐼 0 + 𝛿𝐼. Notably, to study the influence of noise on the stability characteristics, we proceed beyond standard linear stability analysis and consider the Taylor's series expansion up to second order in the perturbation around (𝑆 0 , 𝐼 0 ). To begin with, we rewrite Eq. ( 2) in the following form:

𝜕𝐼 (𝑥, 𝑦, 𝑡) 𝜕𝑡

where

Next, we substitute the perturbations 𝑆 → 𝑆 0 + 𝛿𝑆, and 𝐼 → 𝐼 0 + 𝛿𝐼 into Eq. ( 2) and follow the steps similar to our previous work [13] to obtain the time evolution of moments at an arbitrary grid point in space with discrete correlation function

where 𝐶 𝐼 = 𝜎 2 4Δ𝑥Δ𝑦 for a step size of Δ𝑥 and Δ𝑦 along 𝑥 and 𝑦-axes. Moreover, we simplify the correlators appearing in averaged equations of higher order moments [13] using Novikov's theorem for Gaussian noise process, i.e.,

where

𝛿 𝜉 (𝑠) is the functional derivative which needs to be calculated for non-Markovian noise with correlator defined in Eq. ( 9). Under the assumption that the the correlation time (𝜏) is small, the RHS of Eq. ( 10) can evaluated approximately [27,28]. Following [13], the system of coupled linear equations for the moments may now be written in matrix form:

where 𝑋 = (𝑥 1 , 𝑥 2 , 𝑥 3 , 𝑥 4 , 𝑥 5 ) 𝑇 , and 𝐴 is a 5 × 5 matrix. The components 𝑥 𝑖 are given as

where ℎ = 𝑎 1 + 𝑏 2 , and the other matrix elements are defined as follows:

We note that the equilibrium solution is stable in the absence of diffusion and noise (𝑘 = 0.0, 𝜎 = 0.0). Interestingly, without noise, we only need linear stability to discuss the stability behavior. Therefore, using linear stability, we obtain the condition for Hopf bifurcation as

where (𝑆 0 , 𝐼 0 ) is the equilibrium solution. Without saturation, the critical value of 𝛽 0 at Hopf bifurcation is straightforward to obtain as

ASME Letters in Dynamic Systems and Control / 3

Furthermore, the condition for Turing bifurcation in the absence of noise is obtained from 𝑑𝑒𝑡 ( 𝐴) = 0 and 𝜕 𝜕𝑘 2 [𝑑𝑒𝑡 ( 𝐴)] = 0.

Results

The results are presented under two categories: (1) stability, and

(2) self-organization resulting in pattern formation in the steady state. We obtain stability results by numerically solving for the eigenvalues of 𝐴 in Eq. ( 12) and then plotting the maximal eigenvalue in the plane of the square of the wave number (𝑘 2 ) and the diffusion coefficient (𝐷). This maximal eigenvalue is the leading Lyapunov exponent (local), positive values of which indicate unstable behavior in the neighborhood of the stationary and uniform solutions (𝑆 0 , 𝐼 0 ) [12,13,29]. Furthermore, we obtain the region of instabilities in the parametric plane of the saturation parameter (𝛼) and the correlation time (𝜏) of the O-U noise. In the shaded regions in Figs.

Discussion and Concluding Remarks

We reiterate that higher-order perturbative stability analysis can capture the effects of additive noise that may be latent under lowerorder analyses. Also, contributions from the nonlinear infection force accounting for saturation effects could be more conspicuous under nonlinear stability analysis. Our results indicate that the largest eigenvalue first becomes positive at a certain threshold Fig. 2 The region of instabilities in k 2 -D . The values other of parameters are: b = 1.0, d = 1.0, α = 0.0, β 0 = 35.0, γ = 1.5, C I = 5 × 10 -2 and τ = 0.0.

value of diffusion coefficient, depending upon the saturation and noise parameters. Moreover, the O-U noise correlation time has a notable influence on the spatiotemporal dynamics of the infected population density; precisely, the density shifts to the lower side with an increase in correlation time. Also, notably, the patterns corresponding to lower correlation time approach those obtained for the white noise case (Fig. ( 5)). This offers validation for the numerics since O-U noise, by definition, tends to white noise in the limiting case of vanishing correlation time. Finally, we note that the analysis presented is valid for smaller correlation times. This limitation is related to the evaluation of functional derivatives for a nonlinear system under a correlated noise using Novikov's theorem [27,28]. In light of these results, investigating the emergence of instabilities for other, strongly nonlinear infection forces is a subject of our ongoing research. Furthermore, since nonlinear, stochastic, reaction-diffusion type PDEs arise in a variety of contexts, we expect these results to be significant well beyond epidemic dynamics. Fig. 4 The pattern formation in infected population density in xy plane in the deterministic case. The values of parameters are: b = 1.0, d = 1.0, β 0 = 35.0, γ = 1.5 and D = 10.0. model using COVID-19 spread data [2,3] to be effective in the case of the O-U process as well. Specifically, numerical simulations of spatiotemporal spread -obtained from the PDE epidemic model driven by O-U noise -can be compared against actual infection data in validation studies involving the O-U process. We intend to investigate the aforementioned directions and conclude with the hope that the presented results motivate further research. Fig. 6 The pattern formation in infected population density in xy plane for white noise case with C I = 5 × 10 -2 and τ = 1 × 10 -3 . The values of other parameters are: b = 1.0, d = 1.0, β 0 = 35.0, γ = 1.5 and D = 10.0.