I. INTRODUCTION

Accurate predictive models of epidemic spread -based on real-time infection data -are critical precursors, both for understanding dynamic spread patterns and for enacting interventional control measures for effective mitigation. Despite extensive past research, the COVID-19 pandemic laid bare critical shortcomings in this context and concomitantly highlighted the need to critique both the fundamentals of mathematical modeling of epidemic spread as well as controltheoretic approaches for efficacious interventions. Given that uncertainties -arising, for instance, due to randomness in human behavior and pathogen transmission characteristicscan be critical drivers of epidemic dynamics, the importance of treating epidemics as stochastic dynamic systems is wellrecognized [1]. A fruitful approach in this context is to analyze the spatiotemporal evolution of an epidemic within the framework of a compartmental model represented by a system of coupled partial differential equations (PDE) of the reaction-diffusion type. In previous papers, two of the coauthors (S.R. and M.K.) and their colleagues investigated such a PDE model and validated its predictive capabilities using COVID-19 spread data from the state of Ohio, US [2], [3]. In our more recent work, dynamic instabilities in a related nonlinear PDE model were also uncovered using higher-order perturbation analysis [4]. Here we advance this effort, underscoring the fact that dynamic instabilities that arise both due to nonlinearities, and the interaction between nonlinearity and noise (both of which are justified in the evolution equations of epidemic spread based on physical considerations), can significantly alter the spatiotemporal trajectory of an epidemic as well as the landscape of spread patterns. Furthermore, we note :(1) dynamic instabilities can be potentially correlated with abrupt and drastic changes in epidemic spread patterns (e.g. super spreading events), (2) uncovering instabilities often requires higher-order perturbative analyses, and (3) the study of instabilities and pattern formation in standard reaction-diffusion PDEs, and their profound implications for a broad class of physical systems, has a long history which may be traced back to the remarkable phenomenon known as the Turing instability.

Alan Turing reported that the disparity of diffusion coefficients of the two reacting species involved can drive a homogeneous steady state of a reaction-diffusion PDE system to instability [5], [6], [7]. Experimental evidence for the Turing instability was discovered in the 1990s [8] providing great impetus for further research [5], [9], [10]. A key point is that the classical Turing stability analysis applies only to linear systems. However, motivated by the Turing instability, higher-order perturbative approaches have uncovered the role of nonlinearity in the loss of stability and pattern formation in stochastic systems [9], [11], [12].

We focus on the emergence of instabilities and pattern formation in both the deterministic and stochastic versions of a PDE model of epidemic spread, studied using higherorder perturbative analysis. We now note the significant factors that distinguish the approach and the results in this paper from those in our previously reported work [4]: (1) a different nonlinear function that represents the infection transmission mechanism, taking into account saturation effects (2) instabilities unique to the above nonlinear function that emerge even under an exclusively deterministic analysis (3) new results on noise-induced instabilities, and (4) contour plots of the spatial infection patterns in the steady state, corresponding to numerical solutions of the PDE. We now address the important question of how nonlinearities enter the epidemic PDE model. The PDE model belongs to the class of compartmental epidemic models wherein a population is partitioned into disjoint subsets of Susceptible (S), Infected (I), and Recovered (R) individuals. An epidemic evolves as the I compartment gains individuals due to interaction with the S compartment. Mathematically, the transition between the S and I compartments is represented by the infection force. Supported by empirical evidence of factors such as human behavior patterns and saturation effects, the infection force is often adequately represented by strongly nonlinear functions [13], [14]. As a result, the epidemic model is represented by a system of nonlinear PDEs. We also consider random environmental effects on the infection spread and represent those as a white noise excitation of the Infected population density (I) PDE. We note that the interaction between the nonlinearity and noise in the PDE system can lead to complex dynamics such as stochastic resonance.

The rest of the paper is set as follows. The analytic framework and the research methodology are presented in Section 2. The results are presented in Section 3. The paper concludes in Section 4 with a discussion of the results and concluding remarks.

II. ANALYTIC FRAMEWORK AND METHODOLOGY:

In the complete compartmental model that underlies this work, the full system of reaction-diffusion type PDEs consists of four coupled PDEs -one each for the Susceptible (S), Latent (L), Infected (I), and Recovered (R) population densities corresponding to the four respective compartments [15], [2], [4]. However, since re-infections lie outside the scope of the present analysis, and also since L does not explicitly contribute to the infection force in the model, the dynamics of the reduced system of coupled PDEs for S and I may be meaningfully analyzed for stability. Therefore we focus on this reduced system for our analysis.

We now address the rate of infection. In models without saturation, the rate of infection is defined as β(I) = β 0 I, where β 0 represents the per capita contact, and I is the infected population density. However, taking saturation into account, the rate of infection β(I) can be modeled as [16], [17]

where the term 1 + αI 2 represents the inhibition effect that leads to saturation, i.e. an upper bound on β(I). This is interpreted as a 'psychological' effect reflecting human behavior. This psychological effect is typically a consequence of interventional measures that are enforced -represented by α -such as isolation, quarantine, restriction of public movement, aggressive sanitation, and so on [18]. For lower infection values, a population might take an epidemic less seriously, and this could lead to a rapid increase in the rate of infection. However, as more and more people get infected, individuals are more likely to acknowledge the gravity of the situation and could start responding more positively to protection measures. This behavioral change could be manifest in the increased acceptance of protective measures such as social distancing, sanitation, self-isolation, and masking. This type of behavior indeed influences the rate of further spread and is represented as a nonlinear function, β(I), as specified in Eqn. (1). In this work, α is a nonnegative constant that stands for the rate of saturation. Higher values of α represent faster saturation.

A. Model:

Considering randomness in the infected population density and the nonlinear infection force given in Eqn. (1) [16], [17], the system of coupled PDEs for S and I can be written as:

where b, the birth rate, d, and γ the death rates and D is the diffusion coefficient for the population density S, ∇ 2 is the Laplacian operator, and ξ(t, x, y) is a spatiotemporal Gaussian white noise with

Let (S 0 , I 0 ) be the homogeneous, steady-state solution of the system described by Eq. ( 2). The quantities S 0 , and I 0 are obtained by setting the right-hand side of Eqn. 2 to zero (taking time derivatives to be zero), along with vanishing diffusion and noise terms. Thus (S 0 , I 0 ) simultaneously satisfies:

The acceptable nontrivial solution of Eqn. ( 4) is obtained as

Under no-noise and diffusion-less conditions, this solution (S 0 , I 0 ) is stable, based on sign constraints imposed on the parameter values by physical considerations. However, it is seen that competing populations (i.e. S, and I), if allowed to diffuse, will self-organize in unique patterns. In other words, pattern formation emerges if diffusion causes instability in the homogeneous steady solution -the phenomenon at the heart of the Turing instability. Some points are in order here. Firstly, the ability of the diffusion parameter to create instability could depend on the saturation parameter α. Secondly, it is possible that diffusion alone is insufficient to generate instability. However, diffusing populations may undergo instabilities in the presence of noise of even small intensity. These aspects are investigated in this work.

B. Stability and Moments

For stability analysis, we first perturb the system from its uniform steady-state solution as S 0 → S 0 + δS and I 0 → I 0 + δI. Moreover, in light of our previous results -which underscored the importance of higher-order perturbation analysis -we go beyond standard linear stability analysis and invoke the Taylor series expansion up to second order in the perturbation around (S 0 , I 0 ) in our analysis. Following the analysis similar to our work [4], we write the coupled linear equations for moments in matrix form:

where X = (x 1 , x 2 , x 3 , x 4 , x 5 ) T , and A is a 5 × 5 matrix.

The components x i are given as

where h = a 1 + b 2 , and other matrix elements are defined as:

III. RESULTS

The results fall into two categories -stability, and steadystate pattern formation. The stability results are obtained by numerically solving for the eigenvalues of A in Eqn. (7) and plotting the maximal eigenvalue against the square of the wave number k for each case of interest. We note that this maximal eigenvalue is the leading Lyapunov exponent, positive values for which indicate instability ( [9], [11], [4]). The pattern formation results are obtained by numerically solving for stationary solutions of the PDEs Eqns.2 and obtaining contour plots of the solutions. We note that Eqn. ( 2) is solved using a central difference scheme on a grid size of 100×100, with ∆x = ∆y = 1.0. The time step is ∆t = 0.01, along with no flux boundary conditions. Pattern formation is plotted on the x -y plane in a range of (10, 100) on each axis. The system parameters identical across all cases are b = 1, d = 1.0, γ = 2.8, and β = 35.0. The varying parameters are defined in the corresponding figure captions. In the absence of diffusion D = 0.0 and noise (C I = 0.0), the homogeneous stationary state (S 0 , I 0 ) is stable for the parameter values considered. However, if the susceptible and infected populations are allowed to diffuse, instabilities in the homogeneous stationary state emerge for certain values of the diffusion constant D. Under noise-free conditions, we observe uniquely diffusion-driven instabilities arising for D = 6.0 from the dispersion relation plotted in Fig. (1). The largest eigenvalue (λ) is positive over a finite range of k 2 (square of wave number). This range depends on the value of the saturation parameter (α) and increases with lowering α-values (i.e. lower saturation levels). This diffusion-driven instability produces distinct patterns shown in Figs. ( 4)- (6). As evident from these figures, without saturation (α = 0.0), we obtain a mixed pattern of two types: stripes and dots. However, with increasing values of the saturation parameter, only stripe-like patterns survive [Fig. (6)]. Notably, with a lower value of the diffusion constant D = 3.0, no diffusioninduced instability is observed as the largest eigenvalue is negative for the entire range of k-values [Fig. (2)] (for all the three α values considered); correspondingly, no pattern is seen in Fig. (7). However, if we introduce noise (C I = 0.1) in the infected population dynamics, we see from Fig ( 3) that the largest eigenvalue becomes positive in a finite range of k-values for α = 0.0, and α = 0.5. This is a noisedriven instability and we obtain a noise-induced pattern in Fig. (8) for α = 0.0 with an ensemble size of 10000. While our choice of parameter values is based on numerical experiments, other sets of appropriately chosen values could yield similar results.

IV. DISCUSSION AND CONCLUDING REMARKS

For the nonlinear infection force considered, a first-order perturbation analysis of the PDEs only captured the diffusiondriven instability and was unable to uncover the effects of additive noise on stability. However, a second-order perturbation demonstrates (i) the diffusion driven instabilities (deterministic) and (ii) noise induced instabilities. In the deterministic case, the instability persists for a wider range of k-values, if we lower the value of saturation parameter (α). This range of k-values is smaller for higher levels of saturation. This saturation drives the system from a mixed pattern (stripes and dots) to a only stripes-like pattern, thus, changing the spatial distribution of the infected population in equilibrium. In the stochastic case, noise does not only affect the existing diffusion-driven instability but it also can instigate new instabilities.

To summarize, the results highlight that: (1) nonlinearity and noise, both independently and in tandem, can induce Turing-type instabilities and pattern formation; in particular, both the saturation parameter and the transmission rate are significant factors, (2) higher-order perturbative analysis can uncover such instabilities that might be hidden under lowerorder analyses. Our ongoing research involves the physical characterization of such instabilities with respect to epidemic spread as well as a control-theoretic approach to driving epidemics rapidly to endemic equilibria characterized by specific patterns. Finally, given the wide spectrum of research topics -across disciplines -where reaction-diffusion PDEs and Turing instability play a central role, we expect the results to have broader significance beyond epidemic modeling. Fig. 2. The variation of largest eigenvalue (λ) with the wave number k 2 for second order perturbation. The parameter values are: C I = 0.0, and D = 3.0