1. Introduction
Infectious diseases have always been a focus of global public health security. In recent years, outbreaks of emerging infectious diseases (e.g. HIN1 influenza A) and the drive for population mobility have exacerbated the complexity and uncertainty facing infectious disease prevention and control. The world's low-fertility, aging countries with rising proportions of patients with chronic diseases have exacerbated the healthcare burden and socioeconomic impact. In this context, there is an urgent need for infectious disease research to think from a multidisciplinary cross-cutting perspective: analyzing pathogen adaptations at the biomathematical level, using big data to predict epidemic persistence cycles at the computational mathematical level, and so on. In addition, optimizing public health policies and improving information drivers also provide new ideas for reducing the threat of infectious diseases.
Ordinary differential equations in mathematics, a fundamental tool for dynamical systems [1–3], have deepened their modeling capabilities in the field of infectious disease dynamics through the development of susceptible-infected-recovered (SIR) models [4]. Existing studies have enriched the traditional compartmental framework by introducing mechanisms such as the Allee effect [5], differential mortality weights [6], and nonlinear recovery rates [7]. These analyses extend the traditional framework, while this study combines the core ideas of the literature [6–8]. Other well-recognized features include time-delay effects (e.g. pathogen latency, delayed prevention and control) [9–12]. While our work does not explicitly incorporate time-delay mechanisms, it focuses on improving spatial modeling by constructing a multi-region diffusion–advection framework. In response to the limitations of classical compartmental models in capturing spatio-temporal dynamics, researchers have incorporated spatial variables into epidemic models to enable kinetic analysis [13–15]. For example, Kunz [14] revealed the mechanism of biological pattern formation driven by spatial heterogeneity through coupled diffusion-reaction-convergence, which had broken through the assumption of homogenization under the traditional pattern. Building upon this, network-based epidemic models have been developed to further address the limitations of conventional spatio-temporal frameworks [16–24]. Satorras [16] introduced a networked SIR model based on the Laplace matrix and had revealed the moderating role of the social topology in the propagation dynamics of infectious diseases. Yu [17] further verified the synergistic relationship between network topology and cross-diffusion coefficient on system stability regulation by using the 2SIR model. In another aspect of model development, the advection mechanism has been widely used in various disciplines [25–30]. Wu [29] investigated the long-term dynamics of a cyclic response-advection-dispersal schistosomiasis model, emphasizing the combined effects of advection, seasonality, and spatial heterogeneity on disease transmission. Yang [30] established a network reaction–diffusion model with directional migration term, which revealed the mechanism of pattern formation based on advection effect.
In fact, despite the simplicity of the form of traditional complex network models and their success in capturing the fundamental interactions (pairwise interactions) of complex systems, their fundamental limitation of reducing complex phenomena to a linear superposition of pairwise interactions between nodes has been demonstrated by a large number of studies. Although complex networks effectively model pairwise interactions via edges, numerous socially relevant transmission scenarios (such as shared workplaces and social gatherings) inherently involve synchronous interactions among three or more individuals. In such settings, infection risks stem from the very act of group cohabitation within shared spaces, a characteristic more accurately captured by 2-simplicial complexes than by independent pairwise interactions. The propagation dynamics of real social systems are intrinsically inseparable from higher-order interactions, which has focused the academic community on the construction of higher-order networks [31–38]. Song [31] closed higher-order structures based on dynamic equations break away from traditional node-centered networks, and proposed a higher-order propagation mechanism featuring hyperedge weights. Further theoretical breakthroughs are reflected in the regulation mechanism of the uniformity of the Turing pattern. Through the regulation of the strength of the higher-order interactions in the simple complex form, Luo [35] found that an increase in the strength drived the distribution of the density of infections to homogenization, which was conducive to optimizing the rational allocation of resources during an outbreak of an infectious disease. Building on recent advances in simplicial complex-based modeling—particularly Muolo [36]'s finding that higher-order structures can reshape the conditions for Turing instability—we propose a networked epidemic model that integrates both higher-order diffusion and directional advection. In this framework, higher-order diffusion emerges from three-body interactions defined on a triangular lattice, capturing the coordinated movement among region triplets. Meanwhile, the advection term reflects directional flows influenced by environmental or infrastructural factors. Although these mechanisms have been investigated separately in epidemic dynamics [30, 32], their combined effect on spatiotemporal pattern formation remains largely unexplored. To address this gap and enable scalable simulations, we design a graph-structured dimensionality reduction strategy. This method significantly reduces computational complexity while preserving essential higher-order topological features.
The rest of the paper is organized as follows: In section 2 , we develop a higher-order network-based infectious disease model. In section 3 , we study the necessary conditions for lower-order diffusion, higher-order diffusion and advection to cause Turing instability. In section 4 , we focus on the application of our model to real-world data and networks. Using empirical inter-prefectural mobility data for Japan's 47 prefectures, we construct a directed, weighted commuting network and fit the full heterogeneous-parameter PDE system simultaneously to all prefectures. The fitting is performed with physics-informed neural networks (PINNs) [39, 40], enabling the integration of observational data with the governing dynamical equations. We then assess the role of higher-order spatial interactions in this realistic setting by comparing fits with and without higher-order terms, quantifying the changes in RMSE across prefectures. Finally, we present the conclusions of this paper in section 5 .
2. Modeling
In the transmission of infectious diseases, we divide the population into three parts.
| • | •S (Susceptible individual). Target groups for the transmission of infectious diseases that are not infected but are at risk of being infected. |
| • | •I (Infected individual). Groups of people who carry pathogens, contract diseases and spread them. |
| • | •R (Recovered individual). Groups of people who have successfully cleared pathogens, recovered from infectious diseases and gained immunity. |
In the framework of the classical SIR model, there are three states of nodes in space: susceptible individual, infected individual and recovered individual. We assume that the recovered individual is no longer involved in subsequent propagation, i.e. it has no influence on the dynamics of the S − I subsystem. We then construct an networked infectious disease model with higher-order interactions on a complex network space $\left({\rm{\Omega }}=1,\cdots \,,N\right)$ with discrete nodes. The model reflects the change in node population density over time through an ordinary differential equation with the following expression: 1 ). A is the population replenishment rate, d is the natural mortality rate, μ is the disease-induced mortality rate, and γ is the recovery rate. The model describes the following dynamic behavioral or transformational relationships between different populations in an infectious disease transmission setting:
$\begin{eqnarray}\left\{\begin{array}{l}\displaystyle \frac{{\rm{d}}{S}_{i}\left(t\right)}{{\rm{d}}t}=A-{\rm{d}}{S}_{i}\left(t\right)-\displaystyle \frac{{k}_{1}{I}_{i}\left(t\right)+{k}_{2}{I}_{i}^{2}\left(t\right)}{1+\alpha {I}_{i}^{2}\left(t\right)}{S}_{i}\left(t\right)+{\sum }_{p=1}^{P}{\sigma }_{p}{\sum }_{{j}_{1}=1}^{N}\cdots {\sum }_{{j}_{p}=1}^{N}{A}_{i{j}_{1}\cdots {j}_{p}}^{p}{g}_{p}({S}_{{j}_{1}},\cdots \,,{S}_{{j}_{p}},{I}_{{j}_{1}},\cdots \,,{I}_{{j}_{p}}),\\ \displaystyle \frac{{\rm{d}}{I}_{i}\left(t\right)}{{\rm{d}}t}=\displaystyle \frac{{k}_{1}{I}_{i}\left(t\right)+{k}_{2}{I}_{i}^{2}\left(t\right)}{1+\alpha {I}_{i}^{2}\left(t\right)}{S}_{i}\left(t\right)-(d+\mu ){I}_{i}\left(t\right)-\gamma {I}_{i}\left(t\right)+{\sum }_{p=1}^{P}{\sigma }_{p}{\sum }_{{j}_{1}=1}^{N}\cdots {\sum }_{{j}_{p}=1}^{N}{A}_{i{j}_{1}\cdots {j}_{p}}^{p}{g}_{p}({S}_{{j}_{1}},\cdots \,,{S}_{{j}_{p}},{I}_{{j}_{1}},\cdots \,,{I}_{{j}_{p}}),\\ {S}_{i}\left(0\right)\gt 0,\quad {I}_{i}\left(0\right)\gt 0,\quad i\in {\rm{\Omega }},\quad t\in T=(0,\infty ).\end{array}\right.\end{eqnarray}$
We now give the biological interpretation of parameters in model (| ▶Disease transmissibility ($\frac{{k}_{1}{I}_{i}+{k}_{2}{I}_{i}^{2}}{1+\alpha {I}_{i}^{2}}{S}_{i}$). Describe the number of new infections per unit of time in the susceptible individuals. | |
| ▶Growth rate of disease infectivity (k1Ii). Baseline transmission. | |
| ▶Infectivity of the disease (${k}_{2}{I}_{i}^{2}$). Enhanced spread at high infected density. | |
| ▶Suppression effect ($1+\alpha {I}_{i}^{2}$). Saturation effect decreases infection rate with increasing infected individuals. |
In figure 1, we take the equation for the infected individuals as an example. The population state transition network on the left side demonstrates the coupling mechanism within the nonlinear term, where arrows pointing in indicate an increase in that state and pointing out indicates a decrease. The right side visually portrays how the diffusion term acts, i.e. pairwise interactions through 1-simplexes (links), and three-body interactions through 2-simplexes (triangles). In order to ensure the rationality of the model parameters, based on the physical or realistic significance of the actual application scenarios, all parameters are set to be non-negative real numbers.
Figure 1. Schematic diagram of the SIR higher-order infection epidemic model with birth, death, and recovery. |
As described in [36], σp stands for the coupling strength, we assume that the maximum interaction that exists is $\left(P+1\right)$-body interaction, so p belongs to {1, ⋯ , P}. The (p + 1)-body interactions are denoted by the p-th adjacency tensor Ap, which means that (p + 1)-body interactions between nodes i, j1, ⋯ , jp is given by ${A}_{i{j}_{1}\cdots {j}_{p}}^{p}=1$. In order to construct higher-order coupling terms with diffusion-like properties, we assume that there exists a function hp( · ) for each dimension p such that the coupling function ${g}_{p}({S}_{{j}_{1}},\cdots \,,{S}_{{j}_{p}})$ can be written in the form of a state difference: 2 ) becomes
$\begin{eqnarray*}\begin{array}{rcl}{g}_{p}({S}_{{j}_{1}},\cdots \,,{S}_{{j}_{p}},{I}_{{j}_{1}},\cdots \,,{I}_{{j}_{p}}) & = & {h}_{p}({S}_{{j}_{1}},\cdots \,,{S}_{{j}_{p}},{I}_{{j}_{1}},\cdots \,,{I}_{{j}_{p}})\\ & & -{h}_{p}({S}_{i},\cdots {S}_{i},{I}_{i},\cdots {I}_{i}).\end{array}\end{eqnarray*}$
This construction method is not universally applicable. It has, however, been used as a reasonable modeling assumption in studies such as higher-order coupled oscillator networks (see [31]) and reaction–diffusion networks (see [36]). Its advantages include: ensuring that the coupling terms are zero in the steady state, offering a clear physical analogy with diffusion, and facilitating stability analysis and modal expansion. We adopt this assumption here to enable tractable analysis. Mathematically, the definition of gp when xi = xj ensures that the coupling vanishes when all nodes are in the same state $\begin{eqnarray*}\begin{array}{r}{g}_{p}({S}_{{j}_{1}},\cdots \,,{S}_{{j}_{p}},{I}_{{j}_{1}},\cdots \,,{I}_{{j}_{p}})=0.\end{array}\end{eqnarray*}$
In order to focus on higher-order interactions, we do not consider the cross-diffusion term. This simplification allows us to focus only on pairwise (first-order) interactions and three-body (second-order) interactions without loss of generality. At the same time, we define $F({S}_{i},{I}_{i})\triangleq A-{\rm{d}}{S}_{i}-\frac{{k}_{1}{I}_{i}+{k}_{2}{I}_{i}^{2}}{1+\alpha {I}_{i}^{2}}{S}_{i}$ and $G({S}_{i},{I}_{i})\triangleq \frac{{k}_{1}{I}_{i}+{k}_{2}{I}_{i}^{2}}{1+\alpha {I}_{i}^{2}}{S}_{i}-(d+\mu ){I}_{i}-\gamma {I}_{i}$, we have $\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}{S}_{i}}{{\rm{d}}t}=F({S}_{i},{I}_{i})+{\sigma }_{1}{\sum }_{{j}_{1}=1}^{N}{A}_{i{j}_{1}}^{1}\left({h}_{1}^{S}\left({S}_{{j}_{1}}\right)-{h}_{1}^{S}\left({S}_{i}\right)\right)+{\sigma }_{2}{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2}\left({h}_{2}^{S}\left({S}_{{j}_{1}},{S}_{{j}_{2}}\right)-{h}_{2}^{S}\left({S}_{i},{S}_{i}\right)\right),\\ \frac{{\rm{d}}{I}_{i}}{{\rm{d}}t}=G({S}_{i},{I}_{i})+{\sigma }_{1}{\sum }_{{j}_{1}=1}^{N}{A}_{i{j}_{1}}^{1}\left({h}_{1}^{I}\left({I}_{{j}_{1}}\right)-{h}_{1}^{I}\left({I}_{i}\right)\right)+{\sigma }_{2}{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2}\left({h}_{2}^{I}\left({I}_{{j}_{1}},{I}_{{j}_{2}}\right)-{h}_{2}^{I}\left({I}_{i},{I}_{i}\right)\right).\end{array}\right.\quad i\in {\rm{\Omega }}.\end{eqnarray}$
As pointed out in the literature [41], if higher-order interactions are linear, then they are essentially superpositions of pairwise interactions. For ease of analysis, we focus on the quadratic form as a natural extension of linearity. Specifically, all interactions considered in this study are quadratic. The specific functional forms are $\begin{eqnarray*}\begin{array}{r}\begin{array}{l}{h}_{1}^{S}\left({S}_{{j}_{1}}\right)={D}_{1}^{S}{\left({S}_{{j}_{1}}\right)}^{2},\qquad {h}_{1}^{I}\left({I}_{{j}_{1}}\right)={D}_{1}^{I}{\left({I}_{{j}_{1}}\right)}^{2},\\ \\ {h}_{2}^{S}\left({S}_{{j}_{1}},{S}_{{j}_{2}}\right)={D}_{2}^{S}{S}_{{j}_{1}}{S}_{{j}_{2}},\quad {h}_{2}^{I}\left({I}_{{j}_{1}},{I}_{{j}_{2}}\right)={D}_{2}^{I}{I}_{{j}_{1}}{I}_{{j}_{2}}.\end{array}\end{array}\end{eqnarray*}$
After this treatment, model ( $\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}{S}_{i}}{{\rm{d}}t}=F({S}_{i},{I}_{i})+{\sigma }_{1}{D}_{1}^{S}{\sum }_{{j}_{1}=1}^{N}{A}_{i{j}_{1}}^{1}\left({S}_{{j}_{1}}^{2}-{S}_{i}^{2}\right)+{\sigma }_{2}{D}_{2}^{S}{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2}\left({S}_{{j}_{1}}{S}_{{j}_{2}}-{S}_{i}^{2}\right),\\ \frac{{\rm{d}}{I}_{i}}{{\rm{d}}t}=G({S}_{i},{I}_{i})+{\sigma }_{1}{D}_{1}^{I}{\sum }_{{j}_{1}=1}^{N}{A}_{i{j}_{1}}^{1}\left({I}_{{j}_{1}}^{2}-{I}_{i}^{2}\right)+{\sigma }_{2}{D}_{2}^{I}{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2}\left({I}_{{j}_{1}}{I}_{{j}_{2}}-{I}_{i}^{2}\right).\end{array}\right.\quad i\in {\rm{\Omega }}.\end{eqnarray}$
3. Turing instability anaylisis
3.1. Turing instability conditions for model (3)
This study starts from the stability analysis of an isolated system, where only local reactions are considered and interactions between units are not taken into account. Under this assumption, the equilibrium point of the model (3 ) must be satisfied $F({S}_{i}^{* },{I}_{i}^{* })=0$ and $G({S}_{i}^{* },{I}_{i}^{* })=0$. We denote the equilibrium point $\left({S}_{1}^{* },{I}_{1}^{* },{S}_{2}^{* },{I}_{2}^{* },\cdots \,,{S}_{N}^{* },{I}_{N}^{* }\right)$ as $\left({S}_{* },{I}_{* },{S}_{* },{I}_{* },\cdots \,,{S}_{* },{I}_{* }\right)$. 4 ) is satisfied, ${S}_{* }=\frac{A-\left(d+\mu +\gamma \right){I}_{* }}{d}\gt 0$ necessarily follows. To analyze the stability of model (3 ), we first compute the Jacobian matrix J0 at the equilibrium point (S*, I*, S*, I*, ⋯ , S*, I*) and subsequently derive its characteristic polynomial. The model (3 ) can be decomposed into N identical 2 subsystems, so that the overall characteristic equation of the full system takes the following form, 3 ) remains stable in the absence of disturbances, all characteristic roots must have negative real parts. A necessary condition for this stability is that ${\rm{tr}}({J}^{0})\lt 0$ and ${\rm{\det }}({J}^{0})\gt 0$. Based on this, the necessary condition for the onset of a Turing bifurcation is given by
$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{I}_{* } & = & \frac{-b+\sqrt{{b}^{2}-4ac}}{2a},\\ {S}_{* } & = & \frac{A-\left(d+\mu +\gamma \right){I}_{* }}{d},\end{array}\end{array}\end{eqnarray*}$
where $a=\left({k}_{2}+d\alpha \right)$ $\left(d+\mu +\gamma \right)$, $b={k}_{1}\left(d+\mu +\gamma \right)\,-A{k}_{2}$, $c=d\left(d+\mu +\gamma \right)-A{k}_{1}$. To ensure the relevance of the solution, S* > 0 and I* > 0, the necessary condition for the existence of a positive equilibrium point, needs to be satisfied. We first give the necessary condition for I* > 0 $\begin{eqnarray}\left\{\begin{array}{l}{k}_{1}\left(d+\mu +\gamma \right)-A{k}_{2}\gt 0,\\ A{k}_{1}-d\left(d+\mu +\gamma \right)\gt 0.\end{array}\right.\end{eqnarray}$
If ${I}_{* }\geqslant \frac{A}{d+\mu +\gamma }$, then $\begin{eqnarray*}\begin{array}{r}\begin{array}{l}\Rightarrow 2A\left({k}_{2}+d\alpha \right)+{k}_{1}\left(d+\mu +\gamma \right)-A{k}_{2}\\ -\sqrt{{\left({k}_{1}\left(d+\mu +\gamma \right)-A{k}_{2}\right)}^{2}-4\left({k}_{2}+d\alpha \right)\left(d+\mu +\gamma \right)\left(d\left(d+\mu +\gamma \right)-A{k}_{1}\right)}\leqslant 0,\\ \Rightarrow {A}^{2}d\alpha +d{\left(d+\mu +\gamma \right)}^{2}\leqslant 0.\end{array}\end{array}\end{eqnarray*}$
This contradicts the fundamental assumption that all parameters remain positive. Consequently, the inequality ${I}_{* }\lt \frac{A}{d+\mu +\gamma }$ must hold universally, implying that whenever equation ( $\begin{eqnarray}{\left({\lambda }_{0}^{2}-{\rm{tr}}\left({J}^{0}\right){\lambda }_{0}+\det \left({J}^{0}\right)\right)}^{N}=0,\end{eqnarray}$
where $\begin{eqnarray*}\begin{array}{rcl}{J}^{0} & = & \left|\begin{array}{cc}{j}_{11}^{0} & {j}_{12}^{0}\\ {j}_{21}^{0} & {j}_{22}^{0}\end{array}\right|\\ & = & \left|\begin{array}{cc}-\displaystyle \frac{A}{{S}_{* }} & \displaystyle \frac{{k}_{1}{S}_{* }-2\left(d+\mu +\gamma \right)}{1+\alpha {I}_{* }^{2}}\\ \displaystyle \frac{A}{{S}_{* }}-d & \displaystyle \frac{-{k}_{1}{S}_{* }+2\left(d+\mu +\gamma \right)}{1+\alpha {I}_{* }^{2}}-\left(d+\mu +\gamma \right)\end{array}\right|,\\ \mathrm{tr}\left({J}^{0}\right) & = & -\displaystyle \frac{\left(1+\alpha {I}_{* }^{2}\right)\left(d+\mu +\gamma +\tfrac{A}{{S}_{* }}\right)+{k}_{1}{S}_{* }-2\left(d+\mu +\gamma \right)}{1+\alpha {I}_{* }^{2}},\\ \det \left({J}^{0}\right) & = & \displaystyle \frac{A\left(1+\alpha {I}_{* }^{2}\right)\left(d+\mu +\gamma \right)+{\rm{d}}{S}_{* }\left({k}_{1}{S}_{* }-2\left(d+\mu +\gamma \right)\right)}{{S}_{* }+\alpha {S}_{* }{I}_{* }^{2}}.\end{array}\end{eqnarray*}$
To ensure that the model ( $\begin{eqnarray}\left\{\begin{array}{l}{k}_{1}\left(d+\mu +\gamma \right)-A{k}_{2}\gt 0,\\ A{k}_{1}-2d\left(d+\mu +\gamma \right)\gt 0,\\ {k}_{1}{S}_{* }-2\left(d+\mu +\gamma \right)\gt 0.\end{array}\right.\end{eqnarray}$
3.2. Turing instability induced by pairwise interaction in model (3)
To determine the necessary conditions for Turing instability induced by pairwise interaction, the parameter σ2 in the model (3 ) is set to zero, and a Taylor expansion is carried out around the equilibrium point (S*, I*, S*, I*, ⋯ , S*, I*). Then we retain the linear term in order to obtain the following linearized system
$\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}\delta {S}_{i}}{{\rm{d}}t}={j}_{11}^{0}\delta {S}_{i}+{j}_{12}^{0}\delta {I}_{i}+2{\sigma }_{1}{D}_{1}^{S}{\sum }_{{j}_{1}=1}^{N}{L}_{i{j}_{1}}^{1}\delta {S}_{{j}_{1}}{S}_{* },\\ \frac{{\rm{d}}\delta {I}_{i}}{{\rm{d}}t}={j}_{21}^{0}\delta {S}_{i}+{j}_{22}^{0}\delta {I}_{i}+2{\sigma }_{1}{D}_{1}^{I}{\sum }_{{j}_{1}=1}^{N}{L}_{i{j}_{1}}^{1}\delta {I}_{{j}_{1}}{I}_{* }.\end{array}\right.\quad i\in {\rm{\Omega }}.\end{eqnarray}$
The first-order Laplacian matrix L1 is defined as ${L}^{1}={\left({L}_{ij}^{1}\right)}_{N\times N}={A}^{1}-D$, where D = diag{k1, ⋯ , kN} is the degree matrix, ki represents the order of node i and A1 is a special adjacency tensor, which we also call the adjacency matrix. For an n-th nonzero eigenvalue Λn, the corresponding set of orthogonal eigenvectors is ${{\rm{\Phi }}}_{n}={\left({{\rm{\Phi }}}_{n}^{1},\cdots \,,{{\rm{\Phi }}}_{n}^{N}\right)}^{{\bf{T}}}$. The eigenvector associated with the zero eigenvalue is the all-ones vector Φ0 = (1, ⋯ , 1)T.To express model (7 ) in a compact form, we introduce a perturbation vector $\zeta ={\left(\delta {S}_{1},\delta {I}_{1},\cdots \,,\delta {S}_{N},\delta {I}_{N}\right)}^{{\bf{T}}}$ and the following equation is obtained that 8 ) onto the eigenspace of L1: 9 ) is expanded in Fourier space, i.e. $\left(\begin{array}{c}\delta {\hat{S}}_{n}\\ \delta {\hat{I}}_{n}\end{array}\right)={\sum }_{k=1}^{N}\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right){{\rm{e}}}^{{\lambda }_{k}^{p}t+{\rm{i}}kx}$, and we have 7 ) can be derived as follows
$\begin{eqnarray}\dot{\zeta }=\left({{\mathbb{I}}}_{N}\otimes {J}^{0}+{\sigma }_{1}{L}^{1}\otimes {J}_{{H}_{1}}\right)\zeta ,\end{eqnarray}$
where ${J}_{{H}_{1}}=\left|\begin{array}{ll}2{S}_{* }{D}_{1}^{S} & 0\\ 0 & 2{I}_{* }{D}_{1}^{I}\end{array}\right|$, and ⨂ denotes the Kronecker product, which has a wide range of applications in higher-order networks. Then we project model ( $\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}\delta {\hat{S}}_{n}}{{\rm{d}}t}={j}_{11}^{0}\delta {\hat{S}}_{n}+{j}_{12}^{0}\delta {\hat{I}}_{n}+2{S}_{* }{\sigma }_{1}{D}_{1}^{S}{{\rm{\Lambda }}}_{n}\delta {\hat{S}}_{n},\\ \frac{{\rm{d}}\delta {\hat{I}}_{n}}{{\rm{d}}t}={j}_{21}^{0}\delta {\hat{S}}_{n}+{j}_{22}^{0}\delta {\hat{I}}_{n}+2{I}_{* }{\sigma }_{1}{D}_{1}^{I}{{\rm{\Lambda }}}_{n}\delta {\hat{I}}_{n},\end{array}\right.\quad n\in {\rm{\Omega }},\end{eqnarray}$
the result of our projection is $\left(\begin{array}{c}\delta {\hat{S}}_{n}\\ \delta {\hat{I}}_{n}\end{array}\right)={\sum }_{i=1}^{N}\left(\begin{array}{c}\delta {S}_{i}{{\rm{\Phi }}}_{n}^{i}\\ \delta {I}_{i}{{\rm{\Phi }}}_{n}^{i}\end{array}\right)$. Subsequently, the terms $\delta {\hat{S}}_{n}$ and $\delta {\hat{I}}_{n}$ in model ($\begin{eqnarray*}\begin{array}{r}{\lambda }_{k}^{p}\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right)=\left({J}^{0}+{\sigma }_{1}{{\rm{\Lambda }}}_{k}{J}_{{H}_{1}}\right)\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right),\quad k\in {\rm{\Omega }}.\end{array}\end{eqnarray*}$
The characteristic equation for model ( $\begin{eqnarray}\displaystyle \prod _{k=1}^{N}\left({\left({\lambda }_{k}^{p}\right)}^{2}-{\rm{tr}}\left({J}_{k}^{1}\right){\lambda }_{k}^{p}+det\left(\displaystyle {J}_{k}^{1}\right)\right)=0,\end{eqnarray}$
where $\begin{eqnarray*}\begin{array}{rcl}{J}_{k}^{1} & = & \left|\begin{array}{cc}{j}_{(1,k)}^{1} & {j}_{(2,k)}^{1}\\ {j}_{(3,k)}^{1} & {j}_{(4,k)}^{1}\end{array}\right|\\ & = & \left|\begin{array}{ll}{j}_{11}^{0}+2{S}_{* }{{\rm{\Lambda }}}_{k}{\sigma }_{1}{D}_{1}^{S} & {j}_{12}^{0}\\ {j}_{21}^{0} & {j}_{22}^{0}+2{I}_{* }{{\rm{\Lambda }}}_{k}{\sigma }_{1}{D}_{1}^{I}\end{array}\right|,\quad k\in {\rm{\Omega }}.\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{c}\begin{array}{rcl}{\rm{tr}}({J}_{k}^{1}) & = & {\rm{tr}}({J}^{0})+{\sigma }_{1}{{\rm{\Lambda }}}_{k}{\rm{tr}}({J}_{{H}_{1}}),\,{B}_{2}=4{S}_{\ast }{I}_{\ast }{\sigma }_{1}^{2}{D}_{1}^{S}{D}_{1}^{I},\\ {B}_{1} & = & 2{\sigma }_{1}\left({S}_{\ast }{D}_{1}^{S}{j}_{22}^{0}+{I}_{\ast }{D}_{1}^{I}{j}_{11}^{0}\right),\\ det({J}_{k}^{1}) & = & {B}_{2}{({{\rm{\Lambda }}}_{k})}^{2}+{B}_{1}{{\rm{\Lambda }}}_{k}+{B}_{0},\,{B}_{0}=det({J}^{0}).\end{array}\end{array}\end{eqnarray*}$
To make sure the existence of an eigenroots with a positive real part, it is necessary that at least one of ${\rm{tr}}({J}_{k}^{1})\lt 0$ and ${\rm{\det }}({J}_{k}^{1})\gt 0$ is not true. Given that ${\rm{tr}}({J}_{{H}_{1}})=2({S}_{* }{D}_{1}^{S}+{I}_{* }{D}_{1}^{I})\gt 0$, we can conclude that ${\rm{tr}}({J}_{k}^{1})\lt {\rm{tr}}({J}^{0})\lt 0$. Thus the necessary conditions for low-order diffusion to induce Turing instability is ${\rm{\det }}({J}_{k}^{1})\lt 0$, i.e. $\begin{eqnarray}\left\{\begin{array}{l}{S}_{* }{D}_{1}^{S}{j}_{22}^{0}+{I}_{* }{D}_{1}^{I}{j}_{11}^{0}\gt 0,\\ {\left({\sigma }_{1}{S}_{* }{D}_{1}^{S}{j}_{22}^{0}+{\sigma }_{1}{I}_{* }{D}_{1}^{I}{j}_{11}^{0}\right)}^{2}-4\det \left({J}^{0}\right)\left({S}_{* }{I}_{* }{\sigma }_{1}^{2}{D}_{1}^{S}{D}_{1}^{I}\right)\gt 0.\end{array}\right.\end{eqnarray}$
3.3. Turing instability induced by higher-order interaction in model (3)
To study the necessary conditions for Turing instability induced by higher-order interaction, we perform a Taylor expansion of model (3 ) around the equilibrium point $\left({S}_{* },{I}_{* },{S}_{* },{I}_{* },\cdots \,,{S}_{* },{I}_{* }\right)$ and retain up to the linear term to obtain a linearized system as follows 12 ) in a compact form by introducing a perturbation vector ζ, 13 ), L1 and ${\widetilde{L}}^{2}$ cannot generally be diagonalized simultaneously. That is, it is not possible to project this 2N × 2N system into the same feature space simultaneously. To address this limitation, we will proceed to analyze a specific coupling topology: the reqular topology.
$\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}\delta {S}_{i}}{{\rm{d}}t}={j}_{11}^{0}\delta {S}_{i}+{j}_{12}^{0}\delta {I}_{i}+2{\sigma }_{1}{D}_{1}^{S}{\sum }_{{j}_{1}=1}^{N}{L}_{i{j}_{1}}^{1}\delta {S}_{{j}_{1}}{S}_{* }+2{\sigma }_{2}{D}_{2}^{S}{\sum }_{{j}_{1}=1}^{N}{\widetilde{L}}_{i{j}_{1}}^{2}\delta {S}_{{j}_{1}}{S}_{* },\\ \frac{{\rm{d}}\delta {I}_{i}}{{\rm{d}}t}={j}_{21}^{0}\delta {S}_{i}+{j}_{22}^{0}\delta {I}_{i}+2{\sigma }_{1}{D}_{1}^{I}{\sum }_{{j}_{1}=1}^{N}{L}_{i{j}_{1}}^{1}\delta {I}_{{j}_{1}}{I}_{* }+2{\sigma }_{2}{D}_{2}^{I}{\sum }_{{j}_{1}=1}^{N}{\widetilde{L}}_{i{j}_{1}}^{2}\delta {I}_{{j}_{1}}{I}_{* }.\end{array}\right.\quad i\in {\rm{\Omega }}.\end{eqnarray}$
As described in [36], the generalized Laplace matrix ${\widetilde{L}}^{2}$ for three-body interactions is defined as follows $\begin{eqnarray*}\begin{array}{r}{\widetilde{L}}_{i{j}_{1}}^{2}=\left\{\begin{array}{ll}-{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2} & i=j,\\ {\sum }_{{j}_{1}=1}^{N}{A}_{ij{j}_{1}}^{2} & i\ne j.\end{array}\right.\end{array}\end{eqnarray*}$
If i ≠ j, ${\widetilde{L}}_{ij}^{2}$ denotes the number of 2-simplices formed by edges i and j. Similar to part B of this section, we write model ( $\begin{eqnarray}\dot{\zeta }=\left({{\mathbb{I}}}_{n}\otimes {J}^{0}+{\sigma }_{1}{L}^{1}\otimes {J}_{{H}_{1}}+{\sigma }_{2}{\widetilde{L}}^{2}\otimes {J}_{{H}_{2}}\right)\zeta ,\end{eqnarray}$
where ${J}_{{H}_{2}}=\left|\begin{array}{ll}2{S}_{* }{D}_{2}^{S} & 0\\ 0 & 2{I}_{* }{D}_{2}^{I}\end{array}\right|$. We note that in model (In regular topologies with uniform structure, such as triangular lattice networks with periodic boundary conditions, the higher-order Laplacian matrices (e.g. second-order generalized Laplacians) can often be expressed as scaled versions of first-order Laplacians. This algebraic regularity allows model (13 ), originally defined on a 2N × 2N system, to be projected into a common eigenspace and decomposed into N independent 2 × 2 subsystems, greatly simplifying the analytical complexity.
Inspired by the principles of structural regularity and spectral clarity, we adopt a non-weighted, undirected triangular lattice on a toroidal surface as the underlying topology for our model. While simplified, this regular simplicial structure retains key geometric features relevant to higher-order diffusion. Conceptually, it echoes the analytical tractability demonstrated by the Weighted Triangulated Torus introduced by Wang [33], but deliberately omits the weighted and directional components to focus on the core mechanisms of pattern formation under higher-order interactions. This idealized setting enables clearer interpretation of spatial modes while acknowledging the abstraction from real-world complexity.
In figure 2(a), the triangular lattice torus network's diffusion layer is characterized by six principal connection directions for each node (Ri, Cj), which are Vertical Up: (Ri, Cj) → (Ri, Cj−1), Vertical Down: (Ri, Cj) → (Ri, Cj+1), Horizontal Left: (Ri, Cj) → (Ri−1, Cj), Horizontal Right: (Ri, Cj) → (Ri+1, Cj), Diagonal Southwest: (Ri, Cj) → (Ri−1, Cj+1), Diagonal Northeast: (Ri, Cj) → (Ri+1, Cj−1). Boundary nodes follow a periodic offset connection rule to maintain the torus topology. Let us take the first row as an example: Node (R1, C1) connects to (RN, C2) via diagonal northeast, Node (R1, C2) connects to (RN, C3). Other boundary nodes and so on. In figure 2(b), the specific pattern of three-body interactions in the triangular lattice network can be clearly presented. Now, taking the node labeled with a yellow pentagram as an example, only the participation of this node in three-body interactions with neighboring nodes vertically upward (red line), horizontally to the left (green line), and diagonally to the northeast (blue line) is presented.
Figure 2. The triangular lattice torus network. |
For this triangular lattice torus network, each node participates in the formation of six 1-simplices, and each 1-simplice participates in the formation of two 2-simplices. We have Appendix ), 13 ) takes the following form:
$\begin{eqnarray*}\begin{array}{rc}{L}^{1}=\left\{\begin{array}{ll}-6 & i=j\\ 1 & i\ne j\end{array}\right., & {\widetilde{L}}^{2}=\left\{\begin{array}{ll}-12 & i=j\\ 2 & i\ne j\end{array}\right..\end{array}\end{eqnarray*}$
The eigenvalues of the first-order Laplacian matrix of the aforementioned triangular lattice torus network are obtained through mathematical methods ($\begin{eqnarray*}\begin{array}{rcl}{{\rm{\Lambda }}}_{i,j} & = & 2\cos \left(\frac{2\pi i}{{N}_{0}}\right)+2\cos \left(\frac{2\pi j}{{N}_{0}}\right)\\ & & +2\cos \left(\frac{2\pi i}{{N}_{0}}\right)\cos \left(\frac{2\pi j}{{N}_{0}}\right)\\ & & +2\sin \left(\frac{2\pi i}{{N}_{0}}\right)\sin \left(\frac{2\pi j}{{N}_{0}}\right)-6.\end{array}\end{eqnarray*}$
Based on the relationship of ${\widetilde{L}}^{2}=2{L}^{1}$, the eigenvalues of the second-order generalized Laplacian matrix can be obtained. Meanwhile, model ($\begin{eqnarray*}\begin{array}{r}\dot{\zeta }=\left({{\mathbb{I}}}_{n}\displaystyle \otimes {J}^{0}+{L}^{1}\displaystyle \otimes \left({\sigma }_{1}{J}_{{H}_{1}}+2{\sigma }_{2}{J}_{{H}_{2}}\right)\right)\zeta .\end{array}\end{eqnarray*}$
After projecting onto the eigenspace of L1 and expanding the perturbation in the Fourier space, we obtain $\begin{eqnarray*}\begin{array}{r}{\lambda }_{k}^{h}\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right)=\left({J}^{0}+{{\rm{\Lambda }}}_{k}\left({\sigma }_{1}{J}_{{H}_{1}}+2{\sigma }_{2}{J}_{{H}_{2}}\right)\right)\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right),\quad k\in {\rm{\Omega }}.\end{array}\end{eqnarray*}$
Subsequently, the characteristic equation is derived that $\begin{eqnarray}\displaystyle \prod _{k=1}^{N}\left({\left({\lambda }_{k}^{h}\right)}^{2}-{\rm{tr}}\left({J}_{k}^{2}\right){\lambda }_{k}^{h}+\det \left({J}_{k}^{2}\right)\right)=0,\end{eqnarray}$
where $\begin{eqnarray*}\begin{array}{rcl}{J}_{k}^{2} & = & \left|\begin{array}{cc}{j}_{(1,k)}^{2} & {j}_{(2,k)}^{2}\\ {j}_{(3,k)}^{2} & {j}_{(4,k)}^{2}\end{array}\right|\\ & = & \left|\begin{array}{ll}{j}_{11}^{0}+2{S}_{* }{{\rm{\Lambda }}}_{k}\left({\sigma }_{1}{D}_{1}^{S}+2{\sigma }_{2}{D}_{2}^{S}\right) & {j}_{12}^{0}\\ {j}_{21}^{0} & {j}_{22}^{0}+2{I}_{* }{{\rm{\Lambda }}}_{k}\left({\sigma }_{1}{D}_{1}^{I}+2{\sigma }_{2}{D}_{2}^{I}\right)\end{array}\right|,\quad k\in {\rm{\Omega }}.\end{array}\end{eqnarray*}$
$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{\rm{tr}}({J}_{k}^{2}) & = & {\rm{tr}}({J}^{0})+{{\rm{\Lambda }}}_{k}\left({\sigma }_{1}{\rm{tr}}({J}_{{H}_{1}})+{\sigma }_{2}{\rm{tr}}({J}_{{H}_{2}}),,\right)\\ {C}_{2} & = & 4{S}_{* }{I}_{* }\left({\sigma }_{1}{D}_{1}^{S}+2{\sigma }_{2}{D}_{2}^{S}\right)\left({\sigma }_{1}{D}_{1}^{I}+2{\sigma }_{2}{D}_{2}^{I}\right),\\ {C}_{1} & = & 2\left({S}_{* }\left({\sigma }_{1}{D}_{1}^{S}+2{\sigma }_{2}{D}_{2}^{S}\right){j}_{22}^{0}+{I}_{* }\left({\sigma }_{1}{D}_{1}^{I}+2{\sigma }_{2}{D}_{2}^{I}\right){j}_{11}^{0}\right),\\ \det ({J}_{k}^{2}) & = & {C}_{2}{({{\rm{\Lambda }}}_{k})}^{2}+{C}_{1}{{\rm{\Lambda }}}_{k}+{C}_{0},\\ {C}_{0} & = & \det \left({J}^{0}\right).\end{array}\end{array}\end{eqnarray*}$
Given that ${\rm{tr}}({J}_{{H}_{2}})=2({S}_{* }{D}_{2}^{S}+{I}_{* }{D}_{2}^{I})\gt 0$, we can conclude that ${\rm{tr}}({J}_{k}^{2})\lt {\rm{tr}}({J}^{0})\lt 0$. Thus the necessary conditions for higher-order diffusion to induce Turing instability is ${\rm{\det }}({J}_{k}^{2})\lt 0$, i.e. $\begin{eqnarray}\left\{\begin{array}{l}{S}_{* }\left({\sigma }_{1}{D}_{1}^{S}+2{\sigma }_{2}{D}_{2}^{S}\right){j}_{22}^{0}+{I}_{* }\left({\sigma }_{1}{D}_{1}^{I}+2{\sigma }_{2}{D}_{2}^{I}\right){j}_{11}^{0}\gt 0,\\ {\left({S}_{* }\left({\sigma }_{1}{D}_{1}^{S}+2{\sigma }_{2}{D}_{2}^{S}\right){j}_{22}^{0}+{I}_{* }\left({\sigma }_{1}{D}_{1}^{I}+2{\sigma }_{2}{D}_{2}^{I}\right){j}_{11}^{0}\right)}^{2}-4\det \left({J}^{0}\right)\left({S}_{* }{I}_{* }\left({\sigma }_{1}{D}_{1}^{S}+2{\sigma }_{2}{D}_{2}^{S}\right)\left({\sigma }_{1}{D}_{1}^{I}+2{\sigma }_{2}{D}_{2}^{I}\right)\right)\gt 0.\end{array}\right.\end{eqnarray}$
3.4. Turing instability induced by advection mechanisms in model (3)
In figure 3, we show four advection mechanisms. In the structures shown in figures 3(a) and (b), within the stratosphere, each node has an efferent and incoming degree of 1, and each node points directly downward to a node in the same column. Similar to the connectivity pattern in the diffusion layer, the nodes in the last row point to the nodes in the first row of the same column. We define such a network as a triangular lattice torus network with one-way advection. In figures 3(c) and (d), there are two directional migration paths in the stratosphere. The directed migration paths for susceptible individuals and infected individuals may be the same or opposite, taking into account the effect of factors such as policy and information drivers. These two networks are referred to as triangular lattice torus networks of isotropic advection and reverse advection, respectively. Mathematically, the eigenvalues of the upward and downward advection matrices in figure 3(d) are ${\zeta }_{j}^{u}$ and ${\zeta }_{j}^{d}$ (Appendix ), 3 ) takes the following form
$\begin{eqnarray*}\begin{array}{rcl}{\zeta }_{j}^{u} & = & \cos \left(\frac{2\pi j}{{N}_{0}}\right)+{\rm{i}}\sin \left(\frac{2\pi j}{{N}_{0}}\right)-1,\\ {\zeta }_{j}^{d} & = & \cos \left(\frac{2\pi j}{{N}_{0}}\right)-{\rm{i}}\sin \left(\frac{2\pi j}{{N}_{0}}\right)-1.\end{array}\end{eqnarray*}$
With the introduction of advection, model ( $\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}{S}_{i}}{{\rm{d}}t}=F({S}_{i},{I}_{i})+{{\rm{d}}}_{1}{\sum }_{j=1}^{N}{m}_{ij}{S}_{j}+{\sigma }_{1}{D}_{1}^{S}{\sum }_{{j}_{1}=1}^{N}{A}_{i{j}_{1}}^{1}\left({S}_{{j}_{1}}^{2}-{S}_{i}^{2}\right)+{\sigma }_{2}{D}_{2}^{S}{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2}\left({S}_{{j}_{1}}{S}_{{j}_{2}}-{S}_{i}^{2}\right),\\ \frac{{\rm{d}}{I}_{i}}{{\rm{d}}t}=G({S}_{i},{I}_{i})+{{\rm{d}}}_{2}{\sum }_{j=1}^{N}{n}_{ij}{I}_{j}+{\sigma }_{1}{D}_{1}^{I}{\sum }_{{j}_{1}=1}^{N}{A}_{i{j}_{1}}^{1}\left({I}_{{j}_{1}}^{2}-{I}_{i}^{2}\right)+{\sigma }_{2}{D}_{2}^{I}{\sum }_{{j}_{1}=1}^{N}{\sum }_{{j}_{2}=1}^{N}{A}_{i{j}_{1}{j}_{2}}^{2}\left({I}_{{j}_{1}}{I}_{{j}_{2}}-{I}_{i}^{2}\right),\end{array}\right.\quad i\in {\rm{\Omega }}.\end{eqnarray}$
As described in [30], the matrix $M={({m}_{ij})}_{N\times N}$ and the matrix $N={({n}_{ij})}_{N\times N}$ represent the Laplacian matrices of the directed advection networks corresponding to susceptible individuals and infected individuals, respectively. In addition, the parameters d1 and d2 are used to describe the advection strength that satisfies the conditions d1 ≥ 0 and d2 ≥ 0. $\begin{eqnarray*}\begin{array}{rcl}M & = & \left\{\begin{array}{ll}1,\quad & \,\rm{if}\,\,i\ne j\,\,\rm{and\, there\, exists\, a\, directed\, edge\, from}\,j\,\,\rm{to}\,\,i,\\ 0,\quad & \,\rm{if}\,\,i\ne j\,\,\rm{and\, there\, is\, no\, directed\, edge\, from}\,j\,\,\rm{to}\,\,i,\\ -{\sum }_{{j}_{1}\ne j}{m}_{i{j}_{1}},\quad & \,\rm{if}\,\,i=j.\end{array}\right.\\ N & = & \left\{\begin{array}{ll}1,\quad & \,\rm{if}\,\,i\ne j\,\,\rm{and\, there\, exists\, a\, directed\, edge\, from}\,j\,\,\rm{to}\,\,i,\\ 0,\quad & \,\rm{if}\,\,i\ne j\,\,\rm{and\, there\, is\, no\, directed\, edge\, from}\,j\,\,\rm{to}\,\,i,\\ -{\sum }_{{j}_{1}\ne j}{n}_{i{j}_{1}},\quad & \,\rm{if}\,\,i=j.\end{array}\right.\end{array}\end{eqnarray*}$
Figure 3. Advection methods in the triangular lattice torus network. |
Let us perform a Taylor expansion of model (16 ) around the equilibrium point (S*, I*, S*, I*, ⋯ , S*, I*) and only keep the linear terms. We obtain a linearized system as follows
$\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}\delta {S}_{i}}{{\rm{d}}t}={j}_{11}^{0}\delta {S}_{i}+{j}_{12}^{0}\delta {I}_{i}+{{\rm{d}}}_{1}{\sum }_{j=1}^{N}{m}_{ij}\delta {S}_{j}+2{\sigma }_{1}{D}_{1}^{S}{\sum }_{{j}_{1}=1}^{N}{L}_{i{j}_{1}}^{1}\delta {S}_{{j}_{1}}{S}_{* }+2{\sigma }_{2}{D}_{2}^{S}{\sum }_{{j}_{1}=1}^{N}{\widetilde{L}}_{i{j}_{1}}^{2}\delta {S}_{{j}_{1}}{S}_{* },\\ \frac{{\rm{d}}\delta {I}_{i}}{{\rm{d}}t}={j}_{21}^{0}\delta {S}_{i}+{j}_{22}^{0}\delta {I}_{i}+{{\rm{d}}}_{2}{\sum }_{j=1}^{N}{n}_{ij}\delta {I}_{j}+2{\sigma }_{1}{D}_{1}^{I}{\sum }_{{j}_{1}=1}^{N}{L}_{i{j}_{1}}^{1}\delta {I}_{{j}_{1}}{I}_{* }+2{\sigma }_{2}{D}_{2}^{I}{\sum }_{{j}_{1}=1}^{N}{\widetilde{L}}_{i{j}_{1}}^{2}\delta {I}_{{j}_{1}}{I}_{* },\end{array}\right.\quad i\in {\rm{\Omega }}.\end{eqnarray}$
In the following text, we denote ${M}^{S}\triangleq 2{S}_{* }{\sigma }_{1}{D}_{1}^{S}{L}^{1}\,+2{S}_{* }{\sigma }_{2}{D}_{2}^{S}{\widetilde{L}}^{2}+{d}_{1}M$ and ${N}^{I}\triangleq 2{I}_{* }{\sigma }_{1}{D}_{1}^{I}{L}^{1}+2{I}_{* }{\sigma }_{2}{D}_{2}^{I}{\widetilde{L}}^{2}\,+{d}_{2}N$. To facilitate the analysis, we focus on a special case where MS can be diagonalized and MS = cNI, c ≠ 0. Under this condition, the initial 2N × 2N system can be projected into a common eigenspace and further decomposed into N 2 × 2 subsystems, thereby significantly reducing the overall analytical complexity.
We denote the eigenvalues of MS and NI by θ1, ⋯ , θN and ω1, ⋯ , ωN respectively. Furthermore, for all n ∈ Ω, the corresponding eigenvectors Φn of θn and ωn are identical and form an orthonormal basis. Each eigenvalue can be expressed as ${\theta }_{n}={a}_{n}^{S}+{\rm{i}}{b}_{n}^{S}$ and ${\omega }_{n}={a}_{n}^{I}+{\rm{i}}{b}_{n}^{I}$ for any n ∈ Ω. By the properties of the advection matrix, it follows that ${a}_{n}^{S}\leqslant 0$ and ${a}_{n}^{I}\leqslant 0$.
Similar to part C of this section, the system is projected onto the eigenspace of MS, and the perturbation is expanded in the Fourier space, 18 ), the characteristic equation can be obtained that 19 ) admits at least one root with a positive real part, we assume that the mth eigenvalue satisfies this condition. Let ${\lambda }_{(m,1)}^{a}={x}_{1}+{\rm{i}}{y}_{1}$ and ${\lambda }_{(m,2)}^{a}={x}_{2}+{\rm{i}}{y}_{2}$ denote the two corresponding roots. According to the Vieta formula in the complex domain, 20 ),
$\begin{eqnarray}{\lambda }_{k}^{a}\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right)=\left(\begin{array}{cc}{j}_{11}^{0}+{\theta }_{k} & {j}_{12}^{0}\\ {j}_{21}^{0} & {j}_{22}^{0}+{\omega }_{k}\end{array}\right)\left(\begin{array}{c}{c}_{k}^{1}\\ {c}_{k}^{2}\end{array}\right),\quad k\in {\rm{\Omega }}.\end{eqnarray}$
According to equation ( $\begin{eqnarray}\begin{array}{l}\displaystyle \prod _{k=1}^{N}\left({\left({\lambda }_{k}^{a}\right)}^{2}-\left({\rm{tr}}({J}^{0})+{\theta }_{k}+{\omega }_{k}\right){\lambda }_{k}^{a}\right.\\ \left.+\det ({J}^{0})+{j}_{22}^{0}{\theta }_{k}+{j}_{11}^{0}{\omega }_{k}+{\theta }_{k}{\omega }_{k}\right)=0.\end{array}\end{eqnarray}$
To ensure that equation ( $\begin{eqnarray}\begin{array}{rcl}{\lambda }_{(m,1)}^{a}+{\lambda }_{(m,2)}^{a} & = & {\rm{tr}}({J}^{0})+{\theta }_{m}+{\omega }_{m},\\ {\lambda }_{(m,1)}^{a}{\lambda }_{(m,2)}^{a} & = & \left({j}_{11}^{0}+{\theta }_{m}\right)\left({j}_{22}^{0}+{\omega }_{m}\right)-{j}_{12}^{0}{j}_{21}^{0}.\end{array}\end{eqnarray}$
Next, we consider only the real part of equation ($\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}\mathrm{Re}\left({\lambda }_{(m,1)}^{a}+{\lambda }_{(m,2)}^{a}\right) & = & {\rm{tr}}({J}^{0})+{a}_{m}^{S}+{a}_{m}^{I},\\ \mathrm{Re}\left({\lambda }_{(m,1)}^{a}{\lambda }_{(m,2)}^{a}\right) & = & {x}_{1}{x}_{2}-{y}_{1}{y}_{2}=\det ({J}^{0})+{j}_{11}^{0}{a}_{m}^{I}\\ & & +{j}_{22}^{0}{a}_{m}^{S}+\left({a}_{m}^{S}{a}_{m}^{I}-{b}_{m}^{S}{b}_{m}^{I}\right).\end{array}\end{array}\end{eqnarray*}$
Given that ${\rm{tr}}({J}^{0})\lt 0$, ${a}_{m}^{S}\lt 0$ and ${a}_{m}^{I}\lt 0$, we can conclude that $\mathrm{Re}\left({\lambda }_{(m,1)}^{a}+{\lambda }_{(m,2)}^{a}\right)\lt 0$. Thus, the necessary condition for the existence of the positive real part of the eigenvalue is $\begin{eqnarray}\begin{array}{rcl}{\varphi }_{n} & \triangleq & \det ({J}^{0})+{j}_{11}^{0}{a}_{m}^{I}+{j}_{22}^{0}{a}_{m}^{S}+\left({a}_{m}^{S}{a}_{m}^{I}-{b}_{m}^{S}{b}_{m}^{I}\right)\\ & & +{y}_{1}{y}_{2}\lt 0.\end{array}\end{eqnarray}$
This situation differs from the necessary condition $\det ({J}_{k}^{2})\lt 0$ presented in part C of this section.4. Numerical simulation
In this section, we perform numerical simulations of the theoretical part by using the forward Euler method. The network system is constructed by discretizing the time and space dimensions, where the space step Δx = 1 and the time step Δt = 0.01. The reliability of the theoretical part is verified by numerical simulation of the system by using the Neumann boundary on the region T × Ω. The initial condition is set as a small perturbation near the equilibrium state $\left({S}_{* },{I}_{* },{S}_{* },{I}_{* },\cdots \,,{S}_{* },{I}_{* }\right)$, which is manifested as 3 ) can be simplified. Let us first consider the case where i ≠ j ≠ k,
$\begin{eqnarray*}\begin{array}{l}\begin{array}{l}{S}_{i}\left(0,x\right)={S}_{* }\times \left(1+0.0005\times \mathrm{randn}(1)\right),\\ {I}_{i}\left(0,x\right)={I}_{* }\times \left(1+0.0005\times \mathrm{randn}(1)\right),\end{array}\quad i\in {\rm{\Omega }}.\end{array}\end{eqnarray*}$
In numerical simulations of higher-order interactions, it is challenging to construct three-dimensional second-order Laplace matrices L2 when the number of nodes is high. Even with sparse matrix storage, the construction process is still significantly problematic. To this end, we propose a strategy for dimensionality reduction based on graph structure: as shown in figure 2(b), each edge (1-simplexe) ij belongs to two triangles (2-simplexes) ijk1 and ijk2. By this topological property, the double summation term ${\sum }_{j=1}^{N}{\sum }_{k=1}^{N}{A}_{ijk}^{2}\left({S}_{j}{S}_{k}-{S}_{i}^{2}\right)$ in model ($\begin{eqnarray*}\begin{array}{r}\displaystyle \sum _{j=1}^{N}\displaystyle \sum _{k=1}^{N}{A}_{ijk}^{2}\left({S}_{j}{S}_{k}-{S}_{i}^{2}\right)={A}_{ij{k}_{1}}^{2}{S}_{j}{S}_{{k}_{1}}+{A}_{ij{k}_{2}}^{2}{S}_{j}{S}_{{k}_{2}}-12{S}_{i}^{2}.\end{array}\end{eqnarray*}$
Given that the 3-body interactions among nodes i, j, k is characterized by ${A}_{ijk}^{2}=1$, and that ${\widetilde{L}}^{2}=2$ when i ≠ j, we obtain, $\begin{eqnarray*}\begin{array}{l}{A}_{ij{k}_{1}}^{2}{S}_{j}{S}_{{k}_{1}}+{A}_{ij{k}_{2}}^{2}{S}_{j}{S}_{{k}_{2}}-12{S}_{i}^{2}={S}_{j}\left({S}_{{k}_{1}}+{S}_{{k}_{2}}\right)-12{S}_{i}^{2}\\ \quad =\,2{S}_{j}\frac{{S}_{{k}_{1}}+{S}_{{k}_{2}}}{2}-12{S}_{i}^{2}=\displaystyle \sum _{j=1}^{N}{\widetilde{L}}_{ij}^{2}{S}_{j}\frac{\left({S}_{{k}_{1}}+{S}_{{k}_{2}}\right)}{2}-12{S}_{i}^{2}.\end{array}\end{eqnarray*}$
Moreover, when i = j = k, we observe that $\begin{eqnarray*}\begin{array}{r}\begin{array}{l}\displaystyle \sum _{j=1}^{N}\displaystyle \sum _{k=1}^{N}{A}_{ijk}^{2}\left({S}_{j}{S}_{k}-{S}_{i}^{2}\right)=0,\\ \displaystyle \sum _{j=1}^{N}{\widetilde{L}}_{ij}^{2}{S}_{j}\frac{\left({S}_{{k}_{1}}+{S}_{{k}_{2}}\right)}{2}=-12{S}_{i}^{2}.\end{array}\end{array}\end{eqnarray*}$
Thus, for any 1-simplex ij, we have $\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{j=1}^{N}\displaystyle \sum _{k=1}^{N}{A}_{ijk}^{2}\left({S}_{j}{S}_{k}-{S}_{i}^{2}\right) & = & \displaystyle \sum _{j=1}^{N}\displaystyle \sum _{k=1}^{N}{L}_{ijk}^{2}{S}_{j}{S}_{k}\\ & = & \displaystyle \sum _{j=1}^{N}{\widetilde{L}}_{ij}^{2}{S}_{j}\frac{\left({S}_{{k}_{1}}+{S}_{{k}_{2}}\right)}{2}.\end{array}\end{eqnarray}$
In this way we downscale the three-dimensional second-order Laplacian matrix L2 into a two-dimensional second-order generalized Laplacian matrix ${\widetilde{L}}^{2}$. It is worth noting that although in the triangular lattice torus network we have ${\widetilde{L}}^{2}=2{L}^{1}$, the expression forms of pairwise interactions and higher-order interactions are structurally different. The pairwise term typically involves only self-coupling components such as ${\sum }_{j=1}^{N}{L}_{ij}^{1}{\left({S}_{j}\right)}^{2}$, whereas the higher-order term includes cross-node couplings like ${\sum }_{j=1}^{N}{\widetilde{L}}_{ij}^{2}{S}_{j}\frac{\left({S}_{{k}_{1}}+{S}_{{k}_{2}}\right)}{2}$, and so on. These cross terms represent collaborative effects among neighboring nodes, which cannot be captured by simple pairwise diffusion. In particular, terms like ${\sum }_{j=1}^{N}{\widetilde{L}}_{ij}^{2}{S}_{j}\frac{\left({S}_{{k}_{1}}+{S}_{{k}_{2}}\right)}{2}$ in the higher-order interactions contribute to a different mechanism of influence propagation. Therefore, even if the resulting Laplacian matrices share a scaling relationship, the functional roles of these interaction types in the dynamics are distinct, reflecting not just magnitude but also qualitative differences in the coupling structure.In order to highlight the higher-order interactions, only a small number of typical low-order interaction scenarios are selected as controls for the subsequent simulations in this section. All objects plotted in this section are infected individuals. Unless otherwise specified, the color swatches represent the magnitude of the density of infected individuals.
4.1. Influence of higher-order interactions
To investigate the effect of higher-order interactions on Turing instability, we first focus on analyzing their effect on the eigenvalues of the model.
In figure 4, the parameters are fixed as A = 0.8, d = 0.12, μ = 0.02, γ = 0.6, α = 0.01, k1 = 0.07, k2 = 0.9, and ${D}_{1}^{S}={D}_{2}^{S}=1$. Here, λd denotes the maximum real part of the eigenvalues in equations (10 ) and (14 ), i.e. ${\lambda }_{d}=\max \mathrm{Re}({\lambda }_{i})$ (i = 1, 2, functions of Λn), known as the dispersion relation. A positive λd satisfies the necessary condition for Turing instability. The horizontal axis represents the eigenvalues of the first-order Laplacian L1: light blue dots show pairwise interactions, while light red dots include three-body interactions. Higher-order interactions can either promote or suppress Turing patterns. In figure 4(a) (${D}_{1}^{I}=0.2$, ${D}_{2}^{I}=0.1$, σ1 = σ2 = 0.1), they induce instability; in figure 4(b) (${D}_{1}^{I}=0.1$, ${D}_{2}^{I}=1$, σ1 = σ2 = 10), they suppress it. Since Λn is discrete, the results are interpolated into continuous curves for clarity.
Figure 4. Impact of higher-order interactions on λd. |
Next, let us explore the effect of higher-order interactions on Turing patterns formation. Here, we fix the parameters: A = 0.28, d = 0.12, μ = 0.02, γ = 0.072, α = 0.01, k1 = 0.009, k2 = 0.296, ${D}_{1}^{S}={D}_{2}^{S}=1$, ${D}_{1}^{I}={D}_{2}^{I}=0.1$ and σ1 = 1. It can be observed that the density of infected individuals under this set of parameters shows fewer low-density areas (see figures 5–7).
Figure 5. Influence of σ2 on λd and the density evolution of infected individuals. |
Figure 6. Influence of second-order diffusion coefficients on the density evolution of infected individuals at σ2 = 0.0002. |
Figure 7. Impact of higher-order interactions on Turing patterns formation at time 1000. |
In figure 5, we explore the changes in λd and the evolution of infected density as the second-order coupling strength σ2 increases. In figure 5(a). As σ2 increases, the instability band (i.e. the range of Laplacian eigenvalues Λn with λd > 0) becomes narrower, indicating a reduction in the number of unstable modes. And the lower left corner shows a magnified image. In figure 5(b), we observe that as the coupling strength σ2 increases from 0 to 0.0006, the stabilization time of the pattern gradually lengthens.
Immediately thereafter, we will focus on the influence of key parameters in higher-order interactions on the formation of Turing patterns.
In figure 6, we examine the effect of the second-order diffusion coefficients ${D}_{2}^{S}$ and ${D}_{2}^{I}$ on the pattern evolution process. It can be observed that increasing either coefficient leads to a noticeable prolongation of the stabilization time of the spatial pattern. Combined with the results in figure 5(b), this indicates that the presence of higher-order diffusion tends to slow down the convergence to a steady state in the system dynamics.
In figure 7, we define the discrete coefficient $\mathrm{CV}(x)=\sigma /\bar{x}$, where σ and $\bar{x}$ denote the standard deviation and mean, respectively. A small higher-order diffusion term (σ2: 0 → 0.0002) markedly alters the pattern. At t = 1000 (figures 7(a)–(b)), the barred low-density region transforms into a spiky structure with reduced peak height. Figure 7(c) shows infected densities at t = 1000: blue dots (pairwise only) and orange dots (with 3-body interactions) indicate that higher-order interactions lower infection density. Figure 7(d) plots the time evolution of CV(I) from t = 0 to 600, revealing that higher-order interactions reduce spatial heterogeneity and delay convergence, thus affecting both transient dynamics and steady states of infection spread.
4.2. Influence of advection
In the numerical simulations that follow, if not otherwise specified, the advective methods of susceptible individuals and infected individuals are figures 3(a) and (b) for one-way advection, figure 3(c) for isotropic advection, and figure 3(d) for reverse advection. Next we still use the parameters fixed in part A of this section and fix σ2 = 0.0004.
In figure 8, there is a fixed set of parameters. Specifically, A = 0.9, d = 0.24, μ = 0.02, γ = 0.54, α = 0.01, k2 = 0.9, ${D}_{1}^{S}={D}_{2}^{S}=5$, ${D}_{1}^{I}={D}_{I2}=0.1$, σ1 = 1 and σ2 = 0.1. Similar to figure 4, this side λd is defined as the largest real part of the eigenvalues in equations (14 ) and (19 ), i.e. ${\lambda }_{d}=\max \,\rm{Re}\,({\lambda }_{i})$ (where λi is a function of Λn and i = 2, 3). We still choose the eigenvalues of L1 as the horizontal coordinates, and since the value of λ3 in equation (19 ) is also affected by the advection matrix, the image of the dispersion relation of λ3 exhibits an obvious oscillatory pattern. Here, the blue points indicate diffusion relations without advection, while the orange points include advection. It is observed that advection has a dual effect on the formation of Turing patterns. In figure 8(a), we set k1 = 0.25, d1 = 8, d2 = 0. Observations show that the introduction of advection introduces Turing instability. In figure 8(b), we set k1 = 0.22, d1 = 0, d2 = 8 and observe that the introduction of advection suppresses the Turing instability.
Figure 8. Impact of advection on λd. |
In figure 9, we examine how different advection configurations affect φn and infected density. From equation (21 ), advection-induced Turing instability requires φn < 0. The horizontal axis index i denotes spatial perturbation modes (i, j) from the eigenvalues Λi,j of the first-order Laplacian and advection matrices, with j fixed and i varying to represent spatial frequency in one direction. Multiple φn values for each i are shown vertically. In figure 9(a), point colors indicate advection type: blue (none), red (S-only), green (I-only). Solid symbols satisfy φn < 0; hollow do not. Advection on S lowers φn, increasing unstable modes, while advection on I raises φn, reducing instability. Figure 9(b) shows infected density for the same parameters: with sufficient advection, low-density zones shift from point-like to strip-like patterns, indicating that advection modulates both the extent and morphology of spatial heterogeneity.
Figure 9. Investigate the influence of advection methods on the Φn and patterns. |
Next, two typical advective configurations are considered to explore their influence on Turing pattern formation.
Figure 10 shows the infected density evolution over 1200 time units for all 2500 nodes, with a localized zoomed-in view in the upper right corner illustrating the first 100 nodes over 800 time units. Figures 10(a)–(c) apply one-way advection, while figures 10(d)–(f) implement reverse advection. In both cases, the stabilization time of node-level infected density is significantly shortened as the advection intensity increases. Localized observations further indicate that the upward inclination and density of the striped low-density regions tend to increase with stronger advection. These results suggest that increased advection intensity accelerates the transient dynamics and alters the spatial orientation of low-density regions in the infection field.
Figure 10. The evolution process of the infected density under two advection modes with different advection strength. |
Figure 11 presents the scatter plots of infected density at time 1000 and the temporal evolution of CV(I) over 1200 time units. Figures 11(a), (b) correspond to one-way advection for susceptible individuals, and figures 11(c), (d) show reverse advection. In figure 11(a), increasing advection intensity leads to a broader distribution and an enhanced peak of infected density. A similar trend is observed in figure 11(b) under reverse advection. In figure 11(c), CV(I) increases monotonically with advection intensity, accompanied by a shortened stabilization time. Figure 12(d) exhibits a non-monotonic behavior of CV(I) (first increasing and then decreasing), with a similar shortening of the stabilization period.
Figure 11. The evolution of the infected density distribution under two advection modes with different values of σ2. |
Figure 12. Investigate the influence of advection methods on the Turing patterns at time 600. |
These results demonstrate that advection intensity influences both the degree of spatial heterogeneity and the duration of epidemic system. Specifically, stronger advection can lead to more dispersed patterns and faster convergence to steady-state configurations.
Next, let us explore another scenario where we fix the following parameters A = 0.9, d = 0.24, μ = 0.02, γ = 0.54, α = 0.01, k1 = 0.07, k2 = 0.9, ${D}_{1}^{S}={D}_{2}^{S}=3$, ${D}_{1}^{I}={D}_{2}^{I}=0.1$ and σ1 = 1. It can be observed that the density of Infected under this set of parameters shows more low-density areas. We will explore the effects of advection mechanisms in this case.
In figure 12, we analyze the effects of different advective configurations on the formation of the patterns. We observe that, for a certain convective intensity, the infected density pattern evolves from point-like high-density zones to bar-like high-density zones. This suggests that the introduction of advection can expand the spatial extent of high-density regions in the patterns.
We now focus on the influence of advection on the temporal and spatial dynamics of Turing patterns. Figure 13 compares the infected density evolution in systems with and without advection. Figures 13(a), (c) correspond to simulations without advection, while figures 13(b), (d) include reverse advection. Notably, when advection is introduced [figure 13(b)], the peak infected density increases significantly, and the system exhibits mild oscillations even after reaching a quasi-steady state. In contrast to the horizontally ordered distribution in figures 13(c), (d) displays a more irregular structure, where the low-density strips of infected individuals appear tilted to the upper right, forming an asymmetric and spatially complex pattern. This observation suggests that advection may contribute to the emergence of disordered or chaotic structures in the spatial distribution of infected density.
Figure 13. The influence of d2 = 2 on the evolution process of the density for infected individuals. |
4.3. Real data simulation
To verify the validity of the epidemic model developed above, we carry out numerical simulations using weekly influenza A case data collected from medical institutions across Japan's 47 prefectures, spanning from January 1, 2024 to January 16, 2025.
To incorporate regional heterogeneity in transmission dynamics, we consider a heterogeneous PDE system defined on the empirical inter-prefectural commuting network. The higher-order structure of the empirical network is constructed consistently with the generalized definitions introduced in the theoretical framework. In particular, the second-order adjacency tensor A2 is defined by ${A}_{ijk}^{2}=1$, whenever nodes (i, j, k) form a 2-simplex, and ${A}_{ijk}^{2}=0$ otherwise. Based on A2, the generalized second-order Laplacian matrix ${\widetilde{L}}^{2}$ is given by
$\begin{eqnarray*}{\widetilde{L}}_{ij}^{2}=\left\{\begin{array}{ll}-{\sum }_{{j}_{1}=1}^{47}{\sum }_{{j}_{2}=1}^{47}{A}_{i{j}_{1}{j}_{2}}^{2}, & i=j,\\ {\sum }_{{j}_{1}=1}^{47}{A}_{ij{j}_{1}}^{2}, & i\ne j.\end{array}\right.\end{eqnarray*}$
Specifically, for each prefecture i = 1, …, 47, we introduce heterogeneous parameters, i.e.
$\begin{eqnarray*}\begin{array}{l}{A}_{i},\ {d}_{i},\ {\mu }_{i},\ {\gamma }_{i},\ {k}_{(1,i)},\ {k}_{(2,i)},\ {\alpha }_{i},\\ {D}_{(1,i)}^{S},\ {D}_{(2,i)}^{S},\ {D}_{(1,i)}^{I},\ {D}_{(2,i)}^{I},\ {\sigma }_{(1,i)},\ {\sigma }_{(2,i)},\end{array}\end{eqnarray*}$
which represent local demographic and transmission characteristics.All parameters are treated as node-dependent constants and are jointly inferred together with the neural network representation of the solution fields.
For clarity, the heterogeneous system can also be written explicitly for each prefecture i as
$\begin{eqnarray*}\begin{array}{l}\left\{\begin{array}{l}\displaystyle \frac{{\rm{d}}{S}_{i}(t)}{{\rm{d}}t}={A}_{i}-{d}_{i}{S}_{i}(t)-\displaystyle \frac{{k}_{(1,i)}{I}_{i}(t)+{k}_{(2,i)}{I}_{i}^{2}(t)}{1+{\alpha }_{i}{I}_{i}^{2}(t)}\,{S}_{i}(t)\\ \quad +{\sigma }_{(1,i)}{D}_{(1,i)}^{S}{\sum }_{j=1}^{47}{L}_{ij}^{1}{S}_{j}(t)+{\sigma }_{(2,i)}{D}_{(2,i)}^{S}{\sum }_{j=1}^{47}{\tilde{L}}_{ij}^{2}{S}_{j}(t),\\ \displaystyle \frac{{\rm{d}}{I}_{i}(t)}{{\rm{d}}t}=\displaystyle \frac{{k}_{(1,i)}{I}_{i}(t)+{k}_{(2,i)}{I}_{i}^{2}(t)}{1+{\alpha }_{i}{I}_{i}^{2}(t)}\,{S}_{i}(t)-({d}_{i}+{\mu }_{i}+{\gamma }_{i}){I}_{i}(t)\\ \quad +{\sigma }_{(1,i)}{D}_{(1,i)}^{I}{\sum }_{j=1}^{47}{L}_{ij}^{1}{I}_{j}(t)+{\sigma }_{(2,i)}{D}_{(2,i)}^{I}{\sum }_{j=1}^{47}{\tilde{L}}_{ij}^{2}{I}_{j}(t),\end{array}\right.\quad i=1,2,\ldots ,47.\end{array}\end{eqnarray*}$
Figure 14(a) shows a heat distribution based on the normalized number of infected people per medical institution, with the color scale indicating the normalized values. (Data from https://www.naturalearthdata.com/downloads/10m-cultural-vectors/) figure 14(b) depicts the empirical inter-prefectural commuting network, where nodes represent prefectures and directed edges indicate normalized human mobility flows, with the color scale representing connection strength. (Data from https://www.e-stat.go.jp/dbview?sid=0004019526)
Figure 14. Heat map and commuting network for Japan's 47 prefectures. |
While the triangular lattice adopted in the theoretical analysis provides a mathematically elegant and analytically tractable setting, it cannot fully capture the geographical complexity and social heterogeneity of real-world human interactions. Its primary role in our framework is to isolate and visualize the influence of higher-order spatial interactions—particularly directional advection—in a controlled environment, thereby enabling clearer identification of pattern-forming instabilities and spatial modes. To move beyond the limitations of such idealized settings, we extend our model to a more realistic spatial configuration by incorporating empirical inter-prefectural commuting data for Japan's 47 prefectures. In epidemiological modeling, commuting networks are widely recognized as a natural proxy for the pathways along which infectious diseases spread, since regular population flows between regions directly facilitate pathogen transmission. Figure 14(b) visualizes this empirical commuting network, where nodes correspond to prefectures and directed edges represent population flows. Using official statistics, we construct a directed and weighted mobility network in which edge directionality represents net human flows and raw mobility counts are converted to normalized weights via a proportional scaling method [42]:
$\begin{eqnarray*}\begin{array}{r}{w}_{ij}=\frac{{M}_{ij}}{\displaystyle {\sum }_{k}{M}_{ik}},\quad i,j,k\in 1,2,...\,,\,47,\end{array}\end{eqnarray*}$
where Mij denotes the observed commuting volume from prefecture i to j. This normalization ensures that edge weights reflect the proportion of outbound commuters from each prefecture, thereby preserving relative movement intensities. The resulting weighted migration matrix generalizes the unweighted advection matrix used in the theoretical model, embedding it in a geographically grounded and empirically informed structure. This modification enables the model to account for spatial asymmetry, directional spread, and heterogeneous connectivity, thereby bridging the gap between abstract theory and real-world epidemic dynamics.In scientific research and engineering practice, neural networks have been extensively applied to data fitting and the prediction of complex systems [39, 40, 43–47]. We then employ PINNs [39, 40] to fit the full heterogeneous-parameter PDE system to the data from all 47 prefectures simultaneously. This method integrates the observed epidemiological data with the governing partial differential equations, ensuring consistency with both empirical patterns and theoretical dynamics. Unlike traditional curve-fitting methods or purely data-driven neural networks, PINNs impose physical constraints that help avoid overfitting and improve generalizability, especially in sparse or noisy regions of the dataset. We note that, similar to many gradient-based optimisation methods, the training process of PINNs may exhibit sensitivity to initial parameter values. However, the primary purpose of employing PINNs here is not to locate a globally optimal parameter set, but rather to fit our proposed higher-order network epidemic model by using commuting networks and real-world epidemic data. Crucially, systems incorporating higher-order interactions demonstrate superior fitting accuracy compared to their pairwise interaction counterparts.
Specifically, the PINN takes as input the concatenation of the node identity vector ${e}_{i}\in {{\mathbb{R}}}^{47}$ (the i-th standard basis vector) and time tn, and outputs the state variables (Si(tn), Ii(tn)) for all prefectures. The first-order diffusion is governed by L1, while the higher-order diffusion is governed by ${\widetilde{L}}^{2}$ constructed from the adjacency tensor A2. For each node i, the diffusion items are defined component-wise as 1 ), which are given explicitly by
$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{{ \mathcal D }}_{S,i} & = & {\sigma }_{(1,i)}{D}_{(1,i)}^{S}\displaystyle \sum _{j=1}^{47}{L}_{ij}^{1}{S}_{j}+{\sigma }_{(2,i)}{D}_{(2,i)}^{S}\displaystyle \sum _{j=1}^{47}{\widetilde{L}}_{ij}^{2}{S}_{j},\\ {{ \mathcal D }}_{I,i} & = & {\sigma }_{(1,i)}{D}_{(1,i)}^{I}\displaystyle \sum _{j=1}^{47}{L}_{ij}^{1}{I}_{j}+{\sigma }_{(2,i)}{D}_{(2,i)}^{I}\displaystyle \sum _{j=1}^{47}{\widetilde{L}}_{ij}^{2}{I}_{j},\end{array}\end{array}\end{eqnarray*}$
where σ(1,i) and σ(2,i) denote the node-dependent first- and second-order coupling strengths, respectively. The reaction items follow the model equations defined in model ($\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{ \mathcal F }({S}_{i},{I}_{i}) & = & {A}_{i}-{d}_{i}{S}_{i}-\frac{({k}_{1,i}{I}_{i}+{k}_{2,i}{I}_{i}^{2}){S}_{i}}{1+{\alpha }_{i}{I}_{i}^{2}},\\ { \mathcal G }({S}_{i},{I}_{i}) & = & \frac{({k}_{1,i}{I}_{i}+{k}_{2,i}{I}_{i}^{2}){S}_{i}}{1+{\alpha }_{i}{I}_{i}^{2}}-({d}_{i}+{\mu }_{i}+{\gamma }_{i}){I}_{i}.\end{array}\end{array}\end{eqnarray*}$
Using automatic differentiation to compute temporal derivatives, the physics residuals at time sample tn are
$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{R}_{1,i}({t}_{n}) & = & {\partial }_{t}{S}_{i}({t}_{n})-({{ \mathcal D }}_{S,i}({t}_{n})+{ \mathcal F }({S}_{i}({t}_{n}),{I}_{i}({t}_{n})),\\ {R}_{2,i}({t}_{n}) & = & {\partial }_{t}{I}_{i}({t}_{n})-({{ \mathcal D }}_{I,i}({t}_{n})+{ \mathcal G }({S}_{i}({t}_{n}),{I}_{i}({t}_{n})).\end{array}\end{array}\end{eqnarray*}$
The physics-informed loss is defined as the mean squared residual over all nodes and time samples:
$\begin{eqnarray*}{L}_{p}=\frac{1}{47{N}_{t}}\displaystyle \sum _{n=1}^{{N}_{t}}\displaystyle \sum _{i=1}^{47}\left({R}_{1,i}{({t}_{n})}^{2}+{R}_{2,i}{({t}_{n})}^{2}\right),\end{eqnarray*}$
where Nt is the number of temporal collocation points. The data loss is defined as $\begin{eqnarray*}\begin{array}{rcl}{L}_{d} & = & \displaystyle \frac{1}{47{N}_{t}}\displaystyle \sum _{n=1}^{{N}_{t}}\displaystyle \sum _{i=1}^{47}{\left({I}_{i}({t}_{n})-{I}_{i}^{\mathrm{obs}}({t}_{n})\right)}^{2}\\ & & +\varepsilon \displaystyle \frac{1}{47{N}_{t}}\displaystyle \sum _{n=1}^{{N}_{t}}\displaystyle \sum _{i=1}^{47}{S}_{i}{({t}_{n})}^{2},\end{array}\end{eqnarray*}$
where ϵ = 10−3 prevents the unobserved variable from drifting and ${I}_{i}^{{\rm{o}}{\rm{b}}{\rm{s}}}({t}_{n})$ denotes the observed value of Ii(tn) at time tn.The total loss is
$\begin{eqnarray*}{L}_{t}={L}_{d}+{w}_{p}{L}_{p},\end{eqnarray*}$
where wp follows a warm-up schedule.To assess the contribution of higher-order spatial operators to model performance, we compare fits with and without the inclusion of these terms figure 15 illustrates the change in root mean square error (RMSE) across all prefectures; orange bars indicate a reduction in RMSE when higher-order terms are included. Notably, 57.45% of the prefectures experience improved fit accuracy with the inclusion of higher-order interactions, validating their relevance in real-world modeling.
Figure 15. The change in RMSE across all prefectures. |
Figure 16 further presents the weekly distribution of fitting errors across prefectures over 56 weeks in a heatmap format, where darker colours denote lower errors and higher fitting quality. This visualisation clearly demonstrates the significant role of higher-order spatial operators in enhancing spatio-temporal modelling precision.
Figure 16. Heatmap of weekly fit errors across 56 weeks for all 47 prefectures. |
Figure 17 presents a comparison of the fitted curves with actual observed values across all 47 prefectures, with results for typical groups also displayed separately. The high degree of agreement between the curves indicates that the model demonstrates good overall fitting performance.
Figure 17. Fitted curves versus real data for all 47 prefectures. |
In summary, the introduction of higher-order interactions enhances the model's fitting accuracy and consistency across diverse regions.
5. Conclusion
In this paper, we develop a networked reaction–diffusion infectious disease propagation model on a discrete network space, incorporating higher-order interactions and advection mechanisms. The triangular lattice topology used in the theoretical analysis is chosen solely for its analytical tractability and visual clarity, enabling controlled exploration of the interplay between these mechanisms, rather than as a direct representation of real-world contact structures. By analyzing the necessary conditions for Turing instability under three mechanisms, we identify multiple regulatory pathways for transmission dynamics.
On the one hand, three-body interactions significantly prolong the stabilization period of infectious diseases and inhibit the transition of low-density infected zones from punctate to stripe patterns. On the other hand, the advection mechanism exhibits a spatially selective double-edged effect—when high-density infected areas dominate, advection promotes stripe merging to accelerate spread, whereas when low-density areas dominate, it suppresses spread. The synergy of these mechanisms overcomes the limitations of traditional localized models.
Finally, application to a real-world commuting network demonstrates that incorporating higher-order spatial interactions significantly improves the fitting accuracy of the proposed PINN-based model.
The above work deepens our understanding of the mechanisms underlying Turing pattern formation, providing valuable insights for the study of reaction–diffusion systems. Future research will investigate how higher-order network structures influence dynamical processes and renormalization [48, 49], thereby advancing infectious disease prediction models toward high-dimensional spatiotemporal evolution.
Appendix Prooerties of the triangular lattice torus network
Mathematically, Zhu [30] has obtained the eigenvalues of the composite torus network, which we generalize to the triangular lattice torus network and give the eigenvalues of the first-order Laplace matrix. We define a permutation matrix ${{\rm{\Gamma }}}_{{N}_{0}}$ (${{\rm{\Gamma }}}_{{N}_{0}}^{{\bf{T}}}={{\rm{\Gamma }}}_{{N}_{0}}^{-1}$ and ${N}_{0}=\sqrt{N}$), A1 ), we can get six translational transformations A3 ). A2 ), we obtain the eigenvalues of the first-order Laplace matrix as follows
$\begin{eqnarray*}\begin{array}{rcl}{{\rm{\Gamma }}}_{{N}_{0}} & = & {\rm{\Gamma }}\triangleq \left|\begin{array}{cc}{0}_{{N}_{0}-1\times 1} & {E}_{{N}_{0}-1\times 1}\\ {E}_{1} & 0\end{array}\right|\\ & = & \left|\begin{array}{cccccc}0 & 1 & 0 & \cdots \, & 0 & 0\\ 0 & 0 & 1 & \cdots \, & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots \, & 0 & 1\\ 1 & 0 & 0 & \cdots \, & 0 & 0\end{array}\right|.\end{array}\end{eqnarray*}$
The eigenvalues of this matrix are ${\zeta }_{1}^{{\rm{\Gamma }}},\cdots \,,{\zeta }_{{N}_{0}}^{{\rm{\Gamma }}}$ and its corresponding right eigenvectors are ${\beta }_{1},\cdots \,,{\beta }_{{N}_{0}}$ respectively, i.e. ${\rm{\Gamma }}{\beta }_{n}={\zeta }_{n}^{{\rm{\Gamma }}}{\beta }_{n}$ $\left(n=1,\cdots \,,{N}_{0}\right)$. By transposing the matrix, left-multiplying the inverse matrix, etc, we can obtain $\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}^{-1}{\beta }_{n} & = & {\left({\zeta }_{n}^{{\rm{\Gamma }}}\right)}^{-1}{\beta }_{n},\\ {\beta }_{n}^{{\bf{T}}}{{\rm{\Gamma }}}^{-1} & = & {\zeta }_{n}^{{\rm{\Gamma }}}{\beta }_{n}^{{\bf{T}}},\\ {\beta }_{n}^{{\bf{T}}}{\rm{\Gamma }} & = & {\left({\zeta }_{n}^{{\rm{\Gamma }}}\right)}^{-1}{\beta }_{n}^{{\bf{T}}}.\end{array}\end{eqnarray}$
Next we use the Kronecker product to define an operation $\begin{eqnarray*}\begin{array}{rcl}{\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}} & = & {\left({\beta }_{n}^{i}{\beta }_{m}^{j}\right)}_{{N}_{0}\times {N}_{0}}\\ & = & \left|\begin{array}{cccc}{\beta }_{n}^{1}{\beta }_{m}^{1} & {\beta }_{n}^{1}{\beta }_{m}^{2} & \cdots \, & {\beta }_{n}^{1}{\beta }_{m}^{{N}_{0}}\\ {\beta }_{n}^{2}{\beta }_{m}^{1} & {\beta }_{n}^{2}{\beta }_{m}^{2} & \cdots \, & {\beta }_{n}^{2}{\beta }_{m}^{{N}_{0}}\\ \vdots & \vdots & \vdots & \vdots \\ {\beta }_{n}^{{N}_{0}}{\beta }_{m}^{1} & {\beta }_{n}^{{N}_{0}}{\beta }_{m}^{2} & \cdots \, & {\beta }_{n}^{{N}_{0}}{\beta }_{m}^{{N}_{0}}\end{array}\right|.\end{array}\end{eqnarray*}$
Combining with equation ($\begin{eqnarray}\begin{array}{rcl}{\rm{\Gamma }}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right) & = & {\zeta }_{n}^{{\rm{\Gamma }}}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right),\\ {{\rm{\Gamma }}}^{-1}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right) & = & {\left({\zeta }_{n}^{{\rm{\Gamma }}}\right)}^{-1}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right),\\ \left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right){{\rm{\Gamma }}}^{-1} & = & {\zeta }_{m}^{{\rm{\Gamma }}}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right),\\ \left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right){\rm{\Gamma }} & = & {\left({\zeta }_{m}^{{\rm{\Gamma }}}\right)}^{-1}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right),\\ {{\rm{\Gamma }}}^{-1}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right){{\rm{\Gamma }}}^{-1} & = & {\left({\zeta }_{n}^{{\rm{\Gamma }}}\right)}^{-1}{\zeta }_{m}^{{\rm{\Gamma }}}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right),\\ {\rm{\Gamma }}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right){\rm{\Gamma }} & = & {\zeta }_{n}^{{\rm{\Gamma }}}{\left({\zeta }_{m}^{{\rm{\Gamma }}}\right)}^{-1}\left({\beta }_{n}\displaystyle \otimes {\beta }_{m}^{{\bf{T}}}\right).\end{array}\end{eqnarray}$
For any vector α, the linear transformation L1 satisfies $\begin{eqnarray}\begin{array}{rcl}{L}^{1}\circ \alpha & = & {\rm{\Gamma }}\alpha +{{\rm{\Gamma }}}^{-1}\alpha +\alpha {\rm{\Gamma }}+\alpha {{\rm{\Gamma }}}^{-1}+{{\rm{\Gamma }}}^{-1}\alpha {{\rm{\Gamma }}}^{-1}\\ & & +{\rm{\Gamma }}\alpha {\rm{\Gamma }}-6\alpha ,\end{array}\end{eqnarray}$
where ∘ denotes the linear transformation L1 to acting on the vector α. Let us take ${{\rm{\Psi }}}_{i,j}\triangleq {\beta }_{i}\otimes {\beta }_{j}^{{\bf{T}}}$, α = $\Psi$i,j(i, j = 1, 2, ⋯ , N) and bring them into equation ($\begin{eqnarray*}\begin{array}{rcl}{L}^{1}\circ {{\rm{\Psi }}}_{i,j} & = & \left({\zeta }_{i}^{{\rm{\Gamma }}}+{\left({\zeta }_{i}^{{\rm{\Gamma }}}\right)}^{-1}+{\zeta }_{j}^{{\rm{\Gamma }}}+{\left({\zeta }_{j}^{{\rm{\Gamma }}}\right)}^{-1}\right.\\ & & \left.+{\zeta }_{i}^{{\rm{\Gamma }}}{\left({\zeta }_{j}^{{\rm{\Gamma }}}\right)}^{-1}+{\left({\zeta }_{i}^{{\rm{\Gamma }}}\right)}^{-1}{\zeta }_{j}^{{\rm{\Gamma }}}\right){{\rm{\Psi }}}_{i,j}\\ & & -6{{\rm{\Psi }}}_{i,j}\triangleq {{\rm{\Lambda }}}_{i,j}{{\rm{\Psi }}}_{i,j}.\end{array}\end{eqnarray*}$
Combined with equation ($\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{i,j} & = & 2\cos \left(\frac{2\pi i}{{N}_{0}}\right)+2\cos \left(\frac{2\pi j}{{N}_{0}}\right)\\ & & +2\cos \left(\frac{2\pi i}{{N}_{0}}\right)\cos \left(\frac{2\pi j}{{N}_{0}}\right)\\ & & +2\sin \left(\frac{2\pi i}{{N}_{0}}\right)\sin \left(\frac{2\pi j}{{N}_{0}}\right)-6.\end{array}\end{eqnarray}$
Next, we construct the advection matrix and its eigenvalues are derived as an example for a 4 × 4 network. Then the results will be generalized to general N0 × N0 networks.According to permutation matrix Γ, we can derive the advection matrix in figure 18 as follow A2 ), we can obtain the eigenvalues of M and N are
$\begin{eqnarray*}\begin{array}{r}{M}_{4}=\left|\begin{array}{cccc}{{\rm{\Gamma }}}_{4}-{E}_{4} & 0 & 0 & 0\\ 0 & {{\rm{\Gamma }}}_{4}-{E}_{4} & 0 & 0\\ 0 & 0 & {{\rm{\Gamma }}}_{4}-{E}_{4} & 0\\ 0 & 0 & 0 & {{\rm{\Gamma }}}_{4}-{E}_{4}\end{array}\right|.\end{array}\end{eqnarray*}$
Further we can generalize it to $\begin{eqnarray*}\begin{array}{rcl}M & = & \left|\begin{array}{cccc}{{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}} & 0 & \cdots \, & 0\\ 0 & {{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}} & \cdots \, & 0\\ \vdots & \vdots & \cdots \, & \vdots \\ 0 & 0 & \cdots \, & {{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}}\end{array}\right|,\\ N & = & \left|\begin{array}{cccc}{\left({{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}}\right)}^{-1} & 0 & \cdots \, & 0\\ 0 & {\left({{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}}\right)}^{-1} & \cdots \, & 0\\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots \, & {\left({{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}}\right)}^{-1}\end{array}\right|,\end{array}\end{eqnarray*}$
where $\begin{eqnarray*}\begin{array}{r}\begin{array}{l}{{\rm{\Gamma }}}_{{N}_{0}}-{E}_{{N}_{0}}=\left|\begin{array}{cccc}-1 & 1 & \cdots \, & 0\\ 0 & -1 & \cdots \, & 0\\ \vdots & \vdots & \cdots \, & \vdots \\ 1 & 0 & \cdots \, & -1\end{array}\right|.\end{array}\end{array}\end{eqnarray*}$
Combining with equation ($\begin{eqnarray}\left\{\begin{array}{l}{\zeta }_{j}^{{\rm{up}}}={\rm{\cos }}\left(\frac{2\pi j}{{N}_{0}}\right)+{\rm{i}}{\rm{\sin }}\left(\frac{2\pi j}{{N}_{0}}\right)-1,\\ {\zeta }_{j}^{{\rm{down}}}={\rm{\cos }}\left(\frac{2\pi j}{{N}_{0}}\right)-{\rm{i}}{\rm{\sin }}\left(\frac{2\pi j}{{N}_{0}}\right)-1.\end{array}\right.\end{eqnarray}$
Figure 18. 4 × 4 triangular lattice torus network. |