A deep learning framework for predicting the spread of diffusion diseases

1.   Introduction

The early spread dynamics of infectious diseases have never garnered as much attention as they do today. The use of mathematical models to understand the dynamics of epidemiology can be traced back to 1927 when Kermack first introduced the SIR epidemic model ^[1]. In a closed environment, we assume that the total population size $N$ remains constant and employ a set of ordinary differential equations (ODEs) to describe the mutual transitions among three categories of populations. $S(t), I(t), R(t)$ are functions of time $t$, where $S(t)$ represents susceptible individuals, $I(t)$ represents infected individuals, and $R(t)$ represents recovered individuals. The dynamics of these three populations can be expressed as the following equations

The model in Eq (1.1) assumes that births and deaths are not allowed in a closed environment. Parameter $\beta$ represents the disease transmission rate and $\beta S I$ denotes the number of individuals infected with the disease per unit of time. $\alpha$ represents the recovery rate, signifying the number of individuals recovering and transitioning from the infectious category $I(t)$ to the recovered category $R(t)$ per unit of time.

However, diseases do not initially spread uniformly across the region but also diffuse in space over time. Moreover, different diseases have different ways of spreading. For instance, the COVID-19 pandemic, reached 30 countries within 45 days, which saw the virus geographically spreading. Countries with a higher per-capita reported case counts tended to be highly globalized and had a significant degree of population mobility. Additionally, countries with fewer cases are mostly characterized by lower levels of globalization and domestic mobility. This implies that the disease gradually spreads, and the variation in the number of infections changes over time and space ^[2]. Using PDE to study diffuse disease models. Ndolane combines fractional derivative with diffusion term to study the dynamic behavior of fractional SIR model, which provides a new perspective for the application of fractional diffusion term in epidemic models ^[3].

Reaction-diffusion systems are employed to describe the spatial spread of infectious diseases. We assume that the transmission and recovery rates remain constant over time and space, described by the following partial differential equations

where $\Omega$ is a fixed and bounded domain with a smooth boundary $\partial \Omega$ in $\mathbb{R}^{n}$, $\eta$ represents the outward unit normal vector on the boundary. The homogeneous Neumann boundary conditions imply that the considered diffusion SIR system is self-contained and has no migration across the boundary. In reality, populations always involve new births and deaths, so we consider births and deaths within the region $\Omega$. $B$ denotes the number of newborn individuals and $\mu$ represents the natural death rate. Here we assume the death rates are the same for each population category. The constants $d_i (i = 1, 2, 3)$ are diffusion coefficients representing the diffusion capabilities of each population category. $S(x, t), I(x, t), R(x, t)$ are the population quantities dependent on time and space $S(x, t)+I(x, t)+R(x, t) = N$. Since the top two equations in the system do not contain the recovered individuals $R$, we can analyze the dynamics without the third equation in the system.

Understanding the dynamics of disease spread and accurately estimating the transmission and recovery rates are crucial for epidemiological simulation models. However, accurately estimating these rates often require large amounts of expensive data ^[4]. Moreover, estimating parameters in most infectious disease models remains challenging because disease transmission is inherently stochastic and interdependent with other parameters ^[5].

As is well-known, traditional mathematical methods require only limited data to solve problems. These problems involve inferring the transmission dynamics of infectious diseases based on finite disease dynamic data, including the use of finite difference approximations ^[6]. Sindhu et al. discussed the application of classical statistical distributions ^[7] and parametric methods in medical data analysis ^[8,9]. However, traditional numerical methods are expensive and difficult to obtain data, and they especially perform poorly when dealing with noisy data. Therefore, there is an urgent need to explore new algorithms that can consider limited disease data while leveraging computational models to improve the accuracy, reliability, and predictability of disease dynamics and transmission parameters. Shafiq et al. compared the performance of artificial neural networks (ANN) and traditional parametric methods in epidemic data modeling ^[10], while combining maximum likelihood estimation (MLE) with ANN to explore breast cancer survival data analysis ^[11]. The scientific machine learning approach that combines neural networks with physics-informed modeling has emerged as a new tool in situations with limited or even no data. This in turn extends to the study of a wide range of problems in physics, biology, chemistry, and infectious diseases, including fluid dynamics, solid mechanics, and biochemistry, and the ecology of infectious diseases ^[12]. The Physics-Informed Neural Network (PINN) proposed by Raissi et al. proved to be effective in combining data with physical information^[13]. Specifically, Physics-Informed Neural Network (PINN) has been widely used to solve inverse problems or ill-posed problems by incorporating partial differential equations (PDEs) models with sparse data ^[14]. By modeling nonlinear systems with neural networks and including PDEs in the loss function as a soft penalty term, the amount of required data is reduced. This is achieved using existing knowledge (PDE) to constrain the search space of the model. Yazdani et al. successfully embedded systems biology information into neural networks, creating Systems Biology-Informed Neural Networks (SBINN) ^[15]. Using relatively fewer experimental measurements, they inferred the dynamics of unobserved species and unknown model parameters. Raissi et al. initially applied the Physics-Informed Neural Network method to a standard SIR(ODE) infectious disease system to accurately estimate model parameters ^[16]. Subsequently, Shaier et al. extended the application of the Physics-Informed Neural Network method to various variants of SIR models, generating the unified Differential Information Neural Network (DINN) framework ^[17]. Yin et al. used a generalized partial differential approach to embed ODE disease mechanisms in neural networks, seamlessly coupling deep neural networks and dynamic systems and predicting disease states ^[18]. Difonzo et al. verified that even in the presence of challenging kernel functions, Physics-Informed Neural Network (PINN) can effectively solve inverse problems in periodic dynamics ^[19]. Difonzo et al. proposed applying radial basis functions (RBFs) as activation functions in suitably designed Physics-Informed Neural Network (PINN) to solve the inverse problem ^[20]. These frameworks effectively learn about disease transmission and provide estimates for unknown parameters.

Inspired by the Physics-Informed Neural Network (PINN) framework, we propose the RD-DINN method, which is a deep neural network embedded with reaction-diffusion disease information. This model embeds information about the diffusion disease into a neural network to accurately estimate the dynamics of different population categories and unknown parameters. We consider both the direct (forward) and the inverse problems. The direct problem involves inferring unknown dynamics at any time and location related with known initial and boundary conditions, while the inverse problem aims to deduce unknown parameters, dynamics, and model parameters with limited observed data. We analyze benchmark cases of SIR ^[21] and SEIR ^[22] infectious disease models with diffusion terms, demonstrating the robustness and reliability of our method through various network architectures and the addition of different levels of noise perturbations. We also examine the limitations and shortcomings of our proposed approach in terms of accuracy. The structure of this paper is as follows: In Section 2, we review essential background information on infectious diseases and the development of neural networks and Physics-Informed Neural Network (PINN). In Section 3, we introduce the technical approach of our proposed diffusion disease information neural network. Section 4 showcases the performance of our method on direct and inverse problems in two benchmark cases through computational experiments. In Section 5, we verify the feasibility of the model under real data on the new crown outbreak using five cities in Michigan in April 2020 as distributional conditions for the space $x$ ^*. Finally, in Section 6, we summarize and generalize our approach.

^*https://www.michigan.gov/coronavirus/stats/archived-coronavirus-data

2.   Background models and methods

In the field of mathematical biology, PDE has a better understanding and description of infectious disease propagation, epidemiological dynamics, and other related phenomena exhibiting diffusion behavior. However, when solving most coupled nonlinear PDE systems, it becomes quite challenging, and numerical solutions are the preferred approach. Typically, expensive commercial numerical solvers, such as finite difference method (FDM) ^[23] and finite element method (FEM). Moreover, numerical methods have some drawbacks, including high time consumption, repetitiveness, and a lack of autonomy. To better analyze the dynamic behavior of infectious diseases, it is crucial to build models that can quickly predict diseases. Easily, the SIR equation with a diffusion term and an initial function being positive, has only positive for all $t \in (0, T_{max})$. Therefore, there are two scenarios for the disease spread in the spatiotemporal domain: the disease disappears or persists, which leads to disease-free equilibrium and endemic equilibrium one. We also point out the concept of the basic reproduction number $R_{0}$ ^[24]. In the diffusion SIR model, $R_{0}$ is often expressed as $R_0 = \frac{\beta}{\alpha}$ and serves as a crucial indicator measuring the potential of infectious disease transmission. It represents the average number of susceptible individuals per infected person per unit of time. If $R_0$ is greater than 1, an outbreak is likely and otherwise the epidemic may gradually subside. In our two benchmark cases, considering births and deaths, it is expressed as

Theorem 2.1. Let $\beta$ denote the infection rate, $\alpha$ the recovery rate, $B$ the birth rate, and $\mu$ the natural mortality rate, where $\beta$, $\alpha$, $B$, $\mu$ $\in$ $\mathbb{R} > 0$. Then, the basic reproduction number $R_{0}$ of the system at the disease-free equilibrium (DFE) is given by:

Proof. At the disease-free equilibrium, infected individuals $I = 0$. Thus, the system equation simplifies to

At equilibrium, the time derivative is zero and, assuming a uniform spatial distribution, the spatial derivative is also zero. Therefore, the disease-free equilibrium point $S^* = \frac{B}{\mu}$. In the vicinity of the disease-free equilibrium, we can linearize the system. Let $S = S + s$ and $I = i$, where $s$ and $i$ are small perturbations. By substituting these into the infection equation, ignoring the higher-order terms $\beta si$, we get

In the neighborhood of the disease-free equilibrium, the generation rate of infected individuals is $\beta Si$ and the disappearance rate is ($\alpha + \mu) i$. Thus, the basic reproduction number $R_0$ is the ratio of the generation rate to the disappearance rate

The stability analysis of diseases and the estimation of the basic reproduction number are fundamental explanations for disease dynamics, playing a crucial role in predicting disease spread and formulating preventive measures. However, the ability to track computed values of various populations may vary with different parameter combinations, and there could be variations in real data for different demographic groups. Furthermore, all discussed methods require known initial subpopulation sizes.

With the advancement of science and technology, mathematical biology and epidemiology face a significant challenge in developing a coherent framework that enables them to integrate differential equations with the existing data. We introduce deep learning as a powerful approach, adopting a neural network architecture that mimics brain mechanisms ^[25]. The multi-layer neural network structure is the first perceptron layer that initially processes time and space through weighted sums and biases. Subsequently, it undergoes multiple layers of intricate decision-making until reaching the final decision layer. This decision layer generates outputs corresponding to the dynamic compartments $S(x, t), I(x, t)$, and $R(x, t)$. Deep neural networks seek to minimize the regression mean square error relative to weights and biases by incorporating diffusion terms $\frac{\partial}{\partial x^2}$ into the neural network as PDE residuals through automatic differentiation during backpropagation ^[26] using software such as PyTorch, TensorFlow and incorporating gradient descent methods.

2.1.   Finite difference method (FDM)

The finite difference method (FDM) approximates the partial derivatives in partial differential equations using differences, thereby converting a continuous problem into a discrete one. Discretizing space and time enables computing values at grid points via recursive updates.

Theorem 2.2. Consider the spatial interval $x \in \Omega = [0, L]$ divided into $N$ equidistant grid points with step size $\Delta x$, and time discretized into steps $\Delta t$. The spatial and temporal nodes are given by

The time derivative is approximated using forward difference

The spatial second derivative is approximated using central difference

The no-flow boundary conditions $\partial_\eta S = \partial_\eta I = 0$ are discretized as

The initial conditions $S(x, 0) = S_0(x)$ and $I(x, 0) = I_0(x)$ are discretized as

Using the discretized equations, $S_i^{n+1}$ and $I_i^{n+1}$ can be solved iteratively (e.g., explicit or implicit methods). Explicit methods require stability conditions (e.g., CFL condition), while implicit methods are more stable but require solving linear systems.

2.2.   Physics-Informed Neural Network (PINN)

Physics-Informed Neural Network (PINN) is a data-driven algorithm, which was initially introduced as an alternative model for approximating the solution of differential equations and determining model parameters ^[27]. It utilizes any neural network architecture as the primary framework, such as Fully Connected Neural Networks (FNN) and Convolutional Neural Networks (CNN). The activation function and optimization methodology applied in Physics-Informed Neural Network (PINN) are the same as those of the usual deep learning techniques, and the most intriguing part of the algorithm is the loss function, which is unlike the traditional deep learning methods. Physics-Informed Neural Network (PINN) incorporates physical the residuals of differential equations in the loss function to make the output of the neural network to approximate the system of differential equations. Automatic differentiation is employed to perform derivative calculations with respect to input coordinates and discover underlying dynamics system laws.

2.2.1.   Loss function

Here, $u_{\theta}(x_i, t_i)$ is the neural network's prediction at the data point $(x_i, t_i)$, $u_i$ is the true value at that point, and $N$ is the total number of data points. $\mathcal{L}(u; \lambda)$ is an operator expressed by PDE, and $\lambda$ is the unknown parameters. $M$ is the number of points used to calculate PDE residuals. The total loss is the weighted sum of the data fitting loss and the physical information loss.

where $\omega_1$ and $\omega_2$ are hyperparameters that govern the relative importance of data fitting loss and physical information loss in the total loss.

2.2.2.   Parameter inversion

By minimizing the loss function $L_{\text {total }}$, the parameters $\theta$ and $\lambda$ of the neural network are optimized simultaneously.

where $\theta^{*}$ is the optimal parameter of the neural network and $\lambda^{*}$ is the optimal parameter obtained by inversion. We use the Adam optimizer to optimize both $\theta$ and $\lambda$.

2.3.   Residual connections in deep learning

Residual connections ^[28] have become a cornerstone of modern deep neural networks. These connections address the vanishing gradient problem by facilitating the flow of gradients through the network. Figure 1 shows the residual connection adds the input of a layer directly to its output, enabling the training of much deeper architectures. This approach not only enhances gradient propagation but also improves the overall performance and convergence of the model.

2.4.   Reaction-Diffusion Disease Informed Neural Network (RD-DINN)

Inspired by the latest developments in Physics-Informed Deep Learning and the incorporation of hidden physics models, we propose the RD-DINN framework. This framework utilizes concealed disease information to spread information about diffusive diseases. It employs a fully connected neural network with two-dimensional input and multidimensional output to represent surrogate models for the partial differential equation solutions $S(x, t)$, $I(x, t)$, and $R(x, t)$. The RD-DINN loss function incorporates the encoded disease information model, serving as the computational model for infectious diseases introduced in Section 1. Of course, the embedded disease information could also be a variant of the SIR model, which we later validated with the SEIR model.

Drawing inspiration from the neural attention mechanism in computer vision tasks ^[29], we incorporate residual network modules in addition to fully connected layers, which ensure the preservation of the initial data structure and enhance the weight of hidden states to a significant extent while reducing errors. Figure 2 depicts the schematic diagram of the RD-DINN. We use the residual connection network for both forward and inverse problem experiments.

In Figure refRD-DINN, Res-DNN represents a deep neural network with spatial variable $x$ and temporal variable $t$ as inputs, providing disease dynamics and transmission parameters as outputs. 'PDE Loss' denotes the residual loss of the disease model control equation with a diffusion term, $\omega_r$, $\omega_{ic}$, and $\omega_{bc}$ are the weight values of equation residuals, initial and boundary condition residuals in the loss function. $\mathcal{R}(x, t)$ indicates that the output $u(x, t)$ of neural network meets the dynamic behavior of PDE

where $u$ is the prediction solution of neural network at $\{x, t\}$, $\mathcal{L}(u)$ is differential operator of PDE. $\mathcal{I}(x, t)$ and $\mathcal{B}(x, t)$ are the dynamic behavior of neural network prediction solution $u(x, t)$, satisfying initial and boundary conditions

'Data Loss' accounts for the residual loss between the true values corresponding to randomly sampled points in the spatiotemporal domain in the inverse problem and the network predictions. $\omega_{d}$ represents the weight of the data residual in the loss function of the neural network. $\mathcal{U}^{d}(x, t)-\mathcal{U}(x, t)$ represents the difference between the predicted solution and the real solution of the neural network at $\{x, t\}$. We separate the direct and inverse problems as follows:

● For the direct problem, the loss function comprises three parts: PDE residual loss, initial value loss, and boundary loss.

● For the inverse problem, the loss function consists of two parts: The loss between sampled data true values and Res-DNN outputs, and PDE residual loss.

In all benchmark experiments, the Tanh activation function is utilized due to its smooth gradient properties, which effectively address nonlinear terms (e.g., $S*I$) in the equations and accurately capture complex nonlinear dynamics. Both direct and inverse problems are addressed by minimizing the loss function $L$ using the Adam optimizer, which is employed with hyperparameters: Learning rate $\alpha = 1e-03$, momenta $(\beta_{1}, \beta_{2}) = (0.9, 0.999)$, and $\epsilon = 1e-08$. Convergence is achieved when the total loss falls below $1e-04$. This configuration combines adaptive learning rates with momentum optimization to handle the non-convexity of the loss landscape effectively. The neural network fits the data and infers unobserved disease dynamics and transmission parameters by satisfying the PDE system.

3.   RD-DINN resolves diffusion disease model dynamics and identifies equations

3.1.   Direct problem

Noting that parameters can generate the propagation systems of biological diseases, we restrict the parameters within a specific range near their actual values and within an acceptable margin of error. Therefore, we assume a diffusion disease model with transmission rate $\beta$ and recovery rate $\alpha$ obtained from the literature. For the forward problem, we consider known initial conditions $S(x, 0), I(x, 0)$ and Neumann boundary conditions. We aim to obtain the solutions for the unknown compartments $S(x, t), I(x, t)$.

The next crucial step is to constrain the neural network to comply with the dispersed observed values of $S$ and $I$, and the dynamical system of the control equation. It is necessary to incorporate terms corresponding to the observational data and partial differential equation to construct the loss function.

Let $x\in [0, L]$, the density distribution of the unknown compartments $S(x, t)$ and $I(x, t)$ in the spatiotemporal domain satisfies the partial differential equation system

subject to the initial conditions

and Neumann boundary conditions

Here, the Neumann boundary conditions represent the migration of different population groups at the boundary. We assume that the latent solution $u(x, t)$ can be approximated by a deep neural network $u_{\theta}(x, t)$ with parameters $\theta$. The parameters of the neural network $u_{\theta}(x, t)$ can be estimated by minimizing the mean square loss function of the form

where the PDE residual $\mathcal{L}_{f}(\theta)$, the boundary residual $\mathcal{L}_{bc}(\theta)$, and initial residual $\mathcal{L}_{ic}(\theta)$ can be defined as follows:

Here, $\{x_{d}^{i}, t_{d}^{i}\}^{N_{ic, bc}}_{i = 1}$ denotes the random sampled coordinates of the initial and boundary conditions, and $\{x_{f}^{i}, t_{f}^{i}\}^{N}_{i = 1}$ are a set of spatial-temporal points random sampled within the computational domain $[0, L]\; \times\; [0, T]$ for enforcing the partial differential equation. Here, we employ Latin hypercube sampling as our random sampling method to ensure full coverage of the sampling point range ^[30]. The interval $[x_{min}, x_{max}]$ is divided into $N_{ic}$ layers, with one random point per layer selected to define initial values $S\left(x_i, t_0\right), I\left(x_i, t_0\right)$ at $\left\{\left(x_i, t_0\right)\right\}_{i = 1}^{N_{\mathrm{ic}}}$. The interval $\left[t_0, t_{\mathrm{end}}\right]$ is divided into $N_{bc}$ layers, with one random point per layer selected to define boundary values at $\left\{\left(x_{\text {boundary }}, t_j\right)\right\}_{j = 1}^{N_{\mathrm{bc}}}$ for both $x_{min}$ and $x_{max}$.

In the benchmark experiments, reference solutions for the partial differential equations were computed using a high-order finite difference method. For the SIR model, datasets were generated by systematically varying the parameters $\{\beta, \alpha\} \in \{0.1, 0.3, 0.5, 0.7, 0.9\}$, while for the SEIR model, parameter sets $\{\beta, \sigma, \alpha\} \in \{0.6, 0.8, 1.0, 1.2, 1.4\}$ were employed. A comprehensive sensitivity analysis was conducted across these parameter configurations to rigorously evaluate the robustness of the proposed RD-DINN method under diverse conditions.

To demonstrate the reliability of our algorithm, let's focus on two specific benchmark examples. The parameters for these experiments are given in Table 1.

3.1.1.   Direct SIR model

This benchmark experiment is for the direct problem. We start by considering the SIR model with diffusion

subject to the initial conditions

and the Neumann boundary conditions

Here, $\eta$ is the unit normal vector. Boundary condition equal to 0 means the domain with no population migration in and out of the boundary. $d_1$ and $d_2$ represent the diffusion capabilities of susceptible and infected individuals. $\beta$ represents the transmission rate, and $\alpha$ represents the recovery rate. We assume that these coefficients remain constant. Using the RD-DINN method, we predict the latent solution of the equation. The results are given by Figure 3.

We can observe that when $t$ near zero and $x$ near the boundary, the RD-DINN method exhibits a slightly larger deviation compared to the traditional finite difference method, especially at $(x, t)\in [0, 10]\; \times\; 0.25$. It is reasonable to have a slightly larger deviation due to the abrupt change in population at the early stage of the infectious disease outbreak, as predicting the spread of the disease is challenging during the initial phase. RD-DINN produces excellent predictive solutions, with a relative $L_2$ error of 3%.

3.1.2.   Direct SEIR model

To demonstrate the reliability and reproducibility of our proposed method, the second case selects the SEIR infectious disease model with an added incubation period. The exposed population $E(x, t)$ refers to individuals who have been infected with the virus but have not yet developed symptoms or transmitted the virus to other individuals. Generally, researchers consider the latent period and incubation period as the same concept, with an average latent period of 1.9 days ^[31]. Ndolane studied the numerical method of fractional SEIR model, and discussed the potential application of fractional derivative in spatial propagation model. Fractional derivatives can be combined with diffusion term to describe spatial heterogeneity ^[32]. The SEIR model is described by the following equations

where the recovered population $R$ are not involved. The initial conditions and Neumann boundary conditions are given by

where $\eta$ is the unit normal vector. That Neumann boundary conditions equal to 0 indicates no migration of population at the boundary. Parameters $B$, $\mu$ are consistent with those in the SIR model. $\beta$ represents the disease latent rate, $\sigma$ represents the latent exposure rate, and $\alpha$ represents the disease recovery rate.

Figure 4 shows a schematic diagram of SEIR model prediction. Based on these two disease models mentioned above, we can easily observe that the RD-DINN method, leveraging only a minimal number of initial and boundary data points, along with the selective use of limited data coordinates within the domain to enforce the disease control equations, can effectively predict the unknown dynamics of infectious diseases. The initial and boundary conditions of the model are unknown at the beginning in most cases, which requires inverse problem, leveraging only a small portion of data and incompletely known control equations to find suitable model parameters for modeling.

3.1.3.   Sensitivity analysis

Sensitivity analysis is an important tool for assessing the sensitivity of model outputs to changes in input parameters. Small changes in parameters such as infection rate $\beta$, latency rate $\sigma$, and recovery rate $\alpha$ may lead to significant differences in system behavior. $\frac{\partial}{\partial x}$ describes spatial heterogeneity, and $\{\beta, \sigma, \alpha\}$ are the parameters to be analyzed. We perform individual parameter sensitivity analyses by assigning a specific value to each target parameter $p \in \{\beta, \sigma, \alpha\}$ and keeping the other parameters as defaults. The loss $\mathcal{L}$ is minimized using PINN, the MSE between the predicted and true solution on the test set is calculated and the sensitivity is quantified by analyzing the curve of MSE as a function of parameter $p$

where K is the total number of data points in the test set, $\hat{\mathcal{Q}}^{(k)}$ and $\mathcal{Q}^{(k)}$ are the predicted and true observed values of the model for the $k$th sample, respectively.

In Figure 5, the SIR model maintains MSE consistently below 0.003 across parameter ranges, while the SEIR model maintains MSE under 0.01. This demonstrates the robustness of the RD-DINN framework under parameter variations, where its residual-driven mechanism and integration of prior knowledge mitigate sensitivity effects and enhance adaptability to parameter uncertainties.

3.2.   Inverse problem

In this section, we consider the inverse problem to identify the unknown parameters ($\beta$, $\alpha$) contained in the equations and infer each compartment's unknown dynamics. Specifically, given a small amount of measured values $\{S(x_{d}^{i}, t_{d}^{i}), I(x_{d}^{i}, t_{d}^{i})\}^{N}_{i = 1}$, which are dispersed and may contain noise, We wish to infer a solution $S(x, t), I(x, t)$ in the region $\Omega$ that satisfies the diffusion disease model, and determine the unknown parameters of the model. It is worth emphasizing that for any given data point $\{S(x_{d}^{i}, t_{d}^{i}), I(x_{d}^{i}, t_{d}^{i})\}^{N}_{i = 1}$, the disease information conditions that are met during training are unknown.

So what are the parameters $\beta, \alpha$ that best describe the observed data? To answer this question, we apply RD-DINN to infer the unknown parameters $\beta, \alpha$. Importantly, no initial or boundary conditions are required, which is close to many practical applications where obtaining such information is challenging. We conduct the inverse analysis of two diffusion disease models, assuming that the spatiotemporal domain is subject to dispersed and latent noisy measurements. We explore different network architectures and varying noise levels to demonstrate the robustness and reliability of the RD-DINN in terms of accuracy. To achieve this, we employ a deep neural network to approximate the latent solutions $S(x, t), I(x, t)$ and the unknown parameters $\beta, \alpha$. The training parameters for both are shared, and by minimizing the composite loss below, we train a disease information neural network model that can simultaneously determine unknown parameters and unknown dynamics.

Here, $\{x_{d}^{i}, t_{d}^{i}\}^{N}_{i = 1}$ represents a set of spatial and temporal locations randomly sampled in $[0, L]\; \times\; [0, T]$ during each training iteration. The corresponding population values ${S(x_{d}^{i}, t_{d}^{i}), I(x_{d}^{i}, t_{d}^{i})}$ serve as training datas. We enforce these training datas to the structure dictated by the PDE control equations. The number and position of these location points are identical to the training datas.

To validate the effectiveness of the method, we randomly created a training dataset with $N = 3000$ points based on numerical solutions corresponding to $\beta$ and $\alpha$. Subsequently, a 6-layer deep neural network with 50 neurons in each hidden layer was trained using this dataset. The Adam optimizer is utilized to minimize the loss function. After training completion, the network was fine-tuned to predict the potential solutions $S(x, t)$ and $I(x, t)$ and the unknown parameters $\beta, \alpha$. In the following two benchmark cases, we provide an intuitive evaluation of the accuracy of the Diffusion Disease Information Neural Network. Even in the presence of $5\%$, $10\%$ non-correlated noisy data interference, the network still accurately identifies the unknown parameters and dynamics of the diffusion disease model.

To further evaluate the algorithm performance, we systematically studied the total amount of training data, neural network architecture, and different noise concentrations. The main observation here is that the method exhibits strong robustness to the noise level in the data, achieving reasonable identification accuracy even with noise levels as high as $10\%$. This enhanced robustness seems to outperform competing methods such as Gaussian process regression and methods relying on sparse regression ^[33].

3.2.1.   Inverse SIR model

In the inverse problem, two benchmark cases are considered from the forward problem. The difference lies on the fact that the disease control equations, and the initial and boundary conditions are unknown. Only a tiny, scattered, and possibly noisy dataset is given for the hidden states $S(x, t)$ and $I(x, t)$ to infer the unknown parameters and the unknown dynamics of the equation. The form of the model control equation is given by

define $f(x, t)$ and $g(x, t)$ with

Consider the epidemiological parameters $(\beta, \alpha) = (\lambda_1, \lambda_2)$ be unknown, while $d_1, d_2, B$, and $\mu$ are known. The solutions $S(x, t), I(x, t)$ are predicted through the RD-DINN method, and $\lambda_1$ and $\lambda_2$ are shared parameters with the neural network during training. The predicted solutions, reference solutions, and the pointwise absolute error between them are shown in Figure 6. The RD-DINN method accurately learns the latent solutions $S(x, t), I(x, t)$ of the unknown disease system. The $L_{2}$ error for $S$ is $5.74e-04$, and $I$ is $7.89e-03$.

Figure 6 illustrates the dynamics of the DSIR model across the spatiotemporal domain. The basic regeneration number $R_{0} = \frac{\beta B}{\mu (\mu + \alpha)} \approx 0.83 < 1$ under this set of parameters means that the average number of people infected per infected person is less than 1. The spread of the disease will tend to decline. As infected individuals cannot sustain the transmission of the infection, the disease eventually fades away in the population. It can be observed that starting from the disease outbreak, the number of infected individuals $I(x, t)$ gradually approaches 0, signifying the disappearance of the disease by $t = 10$.

When employing the RD-DINN method for predicting disease dynamics, we simultaneously utilize a shared-weight deep neural network $f_{\theta}(x, t), g_{\theta}(x, t)$ to invert the model's two key parameters: The propagation rate $\beta$ and the recovery rate $\alpha$. The relative error is calculated as (3.14).

In Figure 7, we can see that the two parameters fit well, and because of the initialization of the parameters of the neural network, $\beta$ and $\alpha$ start from zero, and as we keep iterating, the parameters are approaching the actual value. Multiple experiments revealed parameter convergence to true values at 20,000 iterations.

We also conduct a systematic study by maintaining the architecture of the neural network unchanged while varying the noise level and the number of sampled data points to demonstrate the excellent robustness of the algorithm. Regarding the noise level in the data, the method is highly robust. We set 0%, 5%, and 10% noise to data points at 1000, 2000, and 3000, and the MSE of the key parameters of the equation increases with the noise concentration, demonstrating that even if the noise pollution is as high as 10%, the corresponding data points can achieve a reasonable recognition accuracy.

We can observe in Table 2 that none of the $L_2$ errors of RD-DINN in parameter identification exceeds 1%, with 10 repeated trials demonstrating consistent error fluctuations within 0.5% and confirming robustness to neural architecture variations. In Table 3, for adding different concentrations of noise data, none of the $L_2$ errors in parameter identification exceeds 2%, as validated through 10 replications with a maximum error deviation of 0.5%, confirming the method's noise resilience.

In Table 4, it is evident that the residual neural network has a remarkable ability to reduce prediction errors. Specifically, when the number of network layers rises, the inversion error of the architectural parameters of the residual neural network shows a marked decrease. This indicates that the residual neural network is highly effective in improving prediction accuracy and maintaining stability in the face of increased network complexity.

3.2.2.   Inverse SEIR model

In this inverse problem, we employ the SEIR model from the second forward problem, where the parameters of the model $\beta, \sigma, \alpha$ and the dynamics are unknown. Only a small and scattered set of observable data with potential noises are provided for the hidden states $S(x, t), E(x, t), I(x, t)$, and $R(x, t)$. Unlike the first inverse problem case, we introduce an inversion for an additional parameter $\sigma$, which implies that we have less knowledge about the governing laws of the disease information presented by the model. Let

Here, consistent with the inverse problem of the SIR model, we denote $f(t, x), g(t, x), h(t, x)$ as surrogate models for the three equations. The disease dynamics parameters $(\beta, \sigma, \alpha) = (\lambda_1, \lambda_2, \lambda_3)$ are all unknown and predicted through a deep neural network for the unknown dynamics of the compartments $S(x, t), E(x, t)$, and $I(x, t)$. During training, $\lambda_1, \lambda_2$, and $\lambda_3$ are shared parameters with the neural network.

For the SEIR model, we give the representation (3.15) for the basic regeneration number, which is similar to the epidemic models thoroughly described by ordinary differential equations.

Figure 8 shows the constant values observed at $t = 20$ in the entire domain, representing their respective steady states. We compute the basic regeneration number $R_{0}\geq 1$ by bringing in the parameters $\beta, \alpha$, a state in which the populations in each partition have reached a stable equilibrium, signaling the persistence of the disease.

In Figure 9, the three-parameter fits are the same as the actual values, and the neural network parameters are initialized with a value of 1. At the beginning of the iteration, the neural network is randomly fitted, and the values of the fitted parameters are in a decreasing trend. The parameters are gradually approaching the actual values as we keep on iterating and they are the same as the actual parameters at the 50,000 epoch.

Table 5 demonstrates that RD-DINN achieves nearly zero $L_{2}$ error in parameter identification across 10 experimental replications, confirming its accuracy under the basic neural network architecture. As shown in Table 6, while the parameter identification error increases with higher noise concentrations, all errors remain below 1% in 10 replicate experiments, demonstrating consistent performance under noisy conditions.

Model predictions exhibit biases during rapid compartment population shifts, particularly in early outbreaks or unforeseen events, reflecting inherent disease unpredictability and estimation challenges.

3.2.3.   Sensitivity analysis

Figure 10 shows that in the SIR inverse problem for $\alpha = 0.5$, the error rises significantly with increasing $\beta$, indicating that the model is sensitive to high propagation rates. $\beta = 0.5$, the error decreases significantly with increasing $\alpha$, indicating that the model is sensitive to low recovery rates. At high $\beta$, disease outbreak dynamics dominate and model predictions are more stable. The SEIR model inversion errors reached an overall order of magnitude of $1e-07$, indicating a low sensitivity to different parameters. Regarding inter-parameter coupling effects (e.g., simultaneous changes in $\beta$ and $\alpha$), this can be extended to multi-dimensional parameter scans in the future.

4.   Method validation

To validate the model in a public health outbreak, we use data on new crown infections in five cities in mid-Michigan (Allegan, Barry, Eaton, Ingham, and Livings-Ton) in April 2020. We view the properties of a consistent and contiguous spatial distribution of city latitudes as a condition on the spatial $x$ in Figure 11, and we determine whether an infected person is a member of the public through daily nucleic acid tests reported to the public to determine whether an individual is infected. We use population density to represent the current status of each category.

We construct the space of disease for the four cities by performing spatial penalization, considering only transmission at latitude (a one-dimensional spatial problem only). Due to the high degree of dispersion in the data, some degree of data smoothing was required; we compare the impact of the smoothing methods on the data, applying a moving average filter using 3, 4, and 5-day windows, and choosing to apply the 4-day window filter twice; smoothing had little impact on the trend of the data and did not negatively affect the data.

We can see in Figure 12 that the city of Livings-Ton shows an increasing trend in the density of susceptible people and a gradual decreasing trend in the density of infected people in April. Since we do not know this city's birth and death rates, the classical PDE-SIR model with a diffusion term is used to validate this set of accurate data.

By viewing the five cities with the same latitude and close to each other as a one-dimensional spatial distribution, we can see in Figure 13 that RD-DINN can reach an error of $5e-03$ in predicting the spatial distribution. The infection rate is $\beta$ of 0.037 with a float of $10\%$ through the RD-DINN inversion parameter. The recovery rate is $0.070$ with a float of $8\%$. By calculating the mostly regeneration number $R_0 = \frac{\beta}{\alpha} < 1$, it can be concluded that the disease receded briefly in April, which is basically in line with the actual data we selected for the state of Michigan. In Figure 14, we verify the model's extrapolation ability, using the first 20 days of real data to train a pre-training model to reason about the growth of various types of population densities in the latter 10 days. It is not difficult to see that the inference stage and the pre-training stage of the RMSE basically maintains the same order of magnitude, which indicates that the model has the ability to reason about the future state of the disease. Here, the inference phase error has some growth because the model also has some empirical error.

5.   Conclusions

In this paper, we introduce a data-driven deep learning approach based on PINN, which is RD-DINN. RD-DINN is a deep neural network method capable of learning various diseases and how they spread, predicting disease progression, and identifying parameters of disease spread. In the forward problem, where the control equation, initial conditions, and boundary conditions are known, they are incorporated into the neural network loss function. The proposed method achieves an absolute error with finite difference methods of $1e-03$. In the inverse problem, where the inversion of model parameters is crucial, and due to the lack of exact solutions for partial differential equations, we extract a random subset of sampling points from numerical solutions as mini-batch data to minimize the equation and mini-batch data losses. The absolute error between the predicted results and numerical solutions reach $1e-03$. Although RD-DINN demonstrates effective predictions in the forward problem, there is still some error compared to high-order finite differences. We are actively working on more advanced methods to encode disease information into deep neural networks in a more advanced way to improve the errors that occur in positive problems. The inverse problem for time-dependent parameters (e.g., a time-varying transmission rate due to interventions or seasonal effects) is also challenging.

This is the first application of the DINN method to diffusion-based infectious disease models, and the results from two sets of benchmark numerical experiments are promising. Traditional methods such as finite difference and finite element have proven effective in deducing model dynamics. However, RD-DINN not only demonstrates satisfactory performance in predicting model dynamics in the forward problem, but also excels in using mini-batch data for the inverse problem, where it successfully inverts the dynamics and unknown parameters of the model. It is essential to acknowledge that, given the lack of exact solutions for the discussed partial differential equations, the mini-batch data are sampled from the numerical solutions obtained through high-order finite differences.

In the two benchmark disease experiments, we investigate the forward problem of predicting the dynamics of diffusion diseases and the inverse problem of estimating unknown parameters. We also assess the robustness of the RD-DINN method under different neural network architectures and noise levels. We validate the model's effectiveness under new crown data from five cities in Michigan. In conclusion, our proposed Reaction-Diffusion Disease Information Neural Network shows promising prospects as a robust and reliable candidate for learning the dynamics and unknown parameters of compartmental models in infectious diseases, particularly in the context of inverse problems.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors are grateful to the National Natural Science Foundation of China (12471505, 12101576).

Conflict of interest

The authors declare there are no conflicts of interest.

Appendix

Figure A1 illustrates the training process and convergence of neural networks in both direct and inverse problems. In the direct problem, the neural network solves for the system state variables using known partial differential equations, initial conditions, and boundary conditions. The total loss and its components, such as the initial condition loss, boundary condition loss, and function loss, gradually decrease as the number of training epochs increases. This trend indicates that the neural network progressively approximates the true solution. In the inverse problem, the neural network infers the governing equations, initial conditions, or boundary conditions from observed data. Similarly, the total loss and its components decrease with increasing training epochs, demonstrating that the neural network not only fits the observed data but also accurately infers unknown system parameters or dynamic behaviors. These results validate the effectiveness of neural networks in both direct and inverse problems, highlighting their powerful capabilities in modeling complex systems and performing parameter inversion. The total loss in a direct problem consists of three components

where $\mathcal{L}_{\mathrm{IC}}$ is mse for initial conditions, $\mathcal{L}_{\mathrm{BC}}$ is mse for boundary conditions, and $\mathcal{L}_{\mathrm{PDE}}$ is mse for PDE residuals. The total loss in a inverse problem consists of two components

where $\mathcal{L}_{\mathrm{Data}}$ is mse for truth and predict data, $\mathcal{L}_{\mathrm{PDE}}$ mse for PDE residuals.