<figure><img src="image_1.png" /></figure>

HAL
open science

# SIR model: Markovian Approach, Finite Difference Method, PINNs (Physics-Informed Neural Networks)

Charbel Zeus Jupiter Mamlankou

<figure><img src="image_2.png" /></figure>

## To cite this version:

Charbel Zeus Jupiter Mamlankou. SIR model: Markovian Approach, Finite Difference Method, PINNs (Physics-Informed Neural Networks). 2024. hal-04855571

HAL Id: hal-04855571

https://hal.science/hal-04855571v1

Preprint submitted on 25 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SIR Model: Markovian Approach, Finite Difference Method, PINNs (Physics-Informed Neural Networks)

Charbel Zeus Jupiter MAMLANKOU<fnref n="1" />

## Abstract

Recent global health challenges have highlighted the critical importance of robust mathematical frameworks for modeling epidemic spread. Among these models, the SIR (Susceptible-Infected-Recovered) system stands as a fundamental reference. This study focuses on a comprehensive comparison of three distinct approaches for solving the SIR model: the Finite Difference Method (FDM), Markov chains, and Physics-Informed Neural Networks (PINNs). Following a detailed modeling phase, each method is employed to solve the SIR differential system while respecting biological and mathematical constraints such as positivity and invariance. The resulting simulations highlight the strengths and limitations of each approach in terms of accuracy, numerical stability, and computational complexity. This comparative analysis provides enriching perspectives for improving epidemic modeling tools while emphasizing the potential of modern techniques, such as PINNs, to integrate complex and heterogeneous dynamics.

**Keywords:** SIR model, Finite Difference Method, PINNs, Markov chains, modeling, epidemic

## Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 Preliminaries</td><td>2</td></tr><tr><td>3 Finite Difference Method (FDM) approximation</td><td>4</td></tr><tr><td>4 Markov chain approximation</td><td>6</td></tr><tr><td>5 Neural Network approximation</td><td>7</td></tr><tr><td>6 Discussion</td><td>10</td></tr><tr><td>7 Conclusion</td><td>11</td></tr></table>

<footnote n="1">¹National School of Mathematical Engineering and Modeling
E-mail: charbelzeusmamlankou@gmail.com</footnote>
# 1. Introduction

Mathematical modeling of epidemics has become a key field in applied mathematics, providing essential tools for understanding, predicting and controlling the spread of infectious diseases. Among the many models proposed, the SIR (Susceptible-Infected-Recovered) model, first introduced by Kermack and McKendrick in 1927 (1), has established itself as a cornerstone for describing the evolution of epidemics within a population. The model was further developed by Anderson and May (1992), who formalized the mathematical theory of infectious disease, incorporating concepts such as basic reproduction rate and herd immunity. The SIR model is based on a framework of three compartments, each representing a category of individuals in the population: susceptible (S), infected (I) and recovered (R). The dynamic interactions between these compartments are described by a system of differential equations. This formulation makes it possible to analyze how an infection spreads in a closed population.

The model was refined by Diekmann and Heesterbeek in 2000, who incorporated more complex epidemiological scenarios, such as births, deaths and reinfections, enabling a better representation of real epidemic dynamics. For example, their work highlighted the importance of the basic reproduction number $R_0$, which determines whether an epidemic can spread in a given population. Recent contributions, notably those of Ferguson et al. (2020) during the COVID-19 pandemic, have demonstrated the continuing relevance and adaptability of compartmental models in modern epidemiology. Their use to simulate different intervention scenarios has been crucial in informing public health policy. The aim of this study is to carry out a comparative analysis of three distinct approaches to simulating the temporal evolution of population compartments in the SIR model for an X disease. By implementing these different methodologies, we aim to assess their respective effectiveness in solving the system of equations that characterizes the dynamics of the SIR model.

# 2. Preliminaries

This section presents the theoretical and methodological tools used to analyze the SIR model, highlighting the fundamental properties of the system and establishing formal results guaranteeing the validity of the solutions.

The dynamics of infectious diseases are modeled (Lapebie) using the SIR system with mass-action transmission, described as:

$$ \begin{cases} \frac{dS}{dt} = -\beta SI, \\ \frac{dI}{dt} = \beta SI - \gamma I, \\ \frac{dR}{dt} = \gamma I, \end{cases} \qquad (2.1) $$

subject to the initial conditions:

$$ S(0) = S_0, \quad I(0) = I_0, \quad R(0) = R_0, $$

and the population constraint:

$$ S(t) + I(t) + R(t) = N, $$

where N is the total population, assumed constant, β the disease transmission rate and γ the mortality rate.
The mathematical underpinning of our analysis is formalized through the following theorems and propositions. These results serve to establish three critical aspects of the problem: the existence and uniqueness of solutions, the preservation of positivity, and the invariance of solutions within a specific set.

**Theorem 2.1** (Existence and Uniqueness of Solutions). The SIR model defined by (2.1) satisfies the conditions for existence and uniqueness of solutions given initial conditions and a total population constraint.

*Proof:* S(t), I(t) and R(t) are of class C¹, so they are locally Lipschitzian on $\mathbb{R}_+^3$. According to the Cauchy-Lipschitz theorem, we deduce the existence and uniqueness of the solution to the problem associated with our system.

**Proposition 2.2** (Positivity and Boundedness). For all $t \ge 0$, the solutions S(t), I(t), R(t) of the SIR model satisfy $0 \le S(t), I(t), R(t) \le N$.

*Proof:* Assume $S(0), I(0), R(0) \ge 0$. For contradiction, suppose S(t) becomes negative at some $t > 0$. Then $\frac{dS}{dt} = -\beta SI$ implies $\frac{dS}{dt} > 0$ as $S < 0$, which contradicts the assumption. Similar reasoning applies to I(t) and R(t).

The total population is conserved: $S(t) + I(t) + R(t) = N$. Since each component is non-negative, $S, I, R \le N$.

**Theorem 2.3** (Invariant Set). The set $\mathcal{D} = \{(S, I, R) \in \mathbb{R}_+^3 | S + I + R = N\}$ is invariant under the dynamics of the SIR model.

*Proof:* The time derivative of $S + I + R$ is:

$$ \frac{d}{dt}(S + I + R) = \frac{dS}{dt} + \frac{dI}{dt} + \frac{dR}{dt}. $$

Substituting the equations:

$$ \frac{d}{dt}(S + I + R) = (-\beta SI) + (\beta SI - \gamma I) + (\gamma I) = 0. $$

Hence, $S + I + R$ remains constant over time.

If $S(0) + I(0) + R(0) = N$, then $S(t) + I(t) + R(t) = N$ for all $t \ge 0$.

We investigate the progression of disease X in region Y with a total population of 2,500 individuals, where 10 individuals were initially infected prior to the study. The disease-specific parameters - the transmission rate ($\beta$) and recovery rate ($\gamma$) - are set to 0.3 and 0.1, respectively. The complete set of initial parameters for our numerical simulations is as follows:

* Total population: $N = 2500$

* Transmission rate: $\beta = 0.3$

* Recovery rate: $\gamma = 0.1$

* Time horizon: $T = 100$ days
* Initial susceptible population: $S_0 = 2490$

* Initial infected population: $I_0 = 10$

In the next section, we present and analyze the simulation results, evaluating their relevance with respect to our initial objectives for individual method. Furthermore, we discuss the implications of our findings, acknowledge the limitations of our approach, and explore potential avenues for future research.

### 3. Finite Difference Method (FDM) approximation

The finite-difference method is a powerful and widely-used numerical approach for solving the differential equations (2.1) governing the SIR model, enabling accurate simulation of the temporal evolution of population compartments while taking complex epidemic dynamics into account.

To approximate solutions, the equations are discretized with a time step $\Delta t$. Denoting the variables at time $t_n = n\Delta t$ as $S^n, I^n, R^n$, the explicit scheme becomes:

$$ 
\begin{aligned}
S^{n+1} &= S^n - \Delta t \beta S^n I^n, \\
I^{n+1} &= I^n + \Delta t (\beta S^n I^n - \gamma I^n), \\
R^{n+1} &= R^n + \Delta t \gamma I^n.
\end{aligned}
$$

The fixed point corresponds to $I = 0$ (end of epidemic). Let's linearize the equations around this fixed point. Assume:

$$ S^n = S^* + \epsilon_S^n, \quad I^n = I^* + \epsilon_I^n, \quad R^n = R^* + \epsilon_R^n, $$

where $S^*, I^*, R^*$ are the fixed-point values, and $\epsilon_S^n, \epsilon_I^n, \epsilon_R^n$ represent the disturbances. For $I^* = 0$, we have:

$$ S^* = S_0, \quad R^* = N - S_0. $$

Substituting in discrete equations and neglecting higher-order terms in $\epsilon$, we get:

$$ 
\begin{aligned}
\epsilon_S^{n+1} &= \epsilon_S^n - \Delta t \cdot \beta S^* \epsilon_I^n, \\
\epsilon_I^{n+1} &= \epsilon_I^n + \Delta t \cdot (\beta S^* \epsilon_I^n - \gamma \epsilon_I^n), \\
\epsilon_R^{n+1} &= \epsilon_R^n + \Delta t \cdot \gamma \epsilon_I^n.
\end{aligned}
$$

The equations for $\epsilon_S^n$ and $\epsilon_I^n$ are coupled, while $\epsilon_R^n$ can be treated separately. We focus on the coupled subsystem $(\epsilon_S, \epsilon_I)$ to analyze stability.

Writing the system in matrix form, we obtain:

$$ \begin{bmatrix} \epsilon_S^{n+1} \\ \epsilon_I^{n+1} \end{bmatrix} = \begin{bmatrix} 1 & -\Delta t \beta S^* \\ 0 & 1 + \Delta t (\beta S^* - \gamma) \end{bmatrix} \begin{bmatrix} \epsilon_S^n \\ \epsilon_I^n \end{bmatrix}. $$

Let's call the transition matrix $A$:

$$ A = \begin{bmatrix} 1 & -\Delta t \beta S^* \\ 0 & 1 + \Delta t (\beta S^* - \gamma) \end{bmatrix}. $$
The stability of the scheme depends on the eigenvalues of A. If all eigenvalues λ satisfy |λ| ≤ 1, the scheme is stable.

The matrix A is upper triangular, so its eigenvalues are its diagonal elements:

$$ \lambda_1 = 1, \quad \lambda_2 = 1 + \Delta t(\beta S^* - \gamma). $$

For |λ₂| ≤ 1, we need :

$$ -1 \le 1 + \Delta t(\beta S^* - \gamma) \le 1. $$

This gives two inequalities:

• 1 + Δt(βS* - γ) ≥ -1, or Δt ≤ 2 / (γ - βS*)

• 1 + Δt(βS* - γ) ≤ 1, or Δt ≥ 0.

The first condition is more restrictive.

Remark 3.1. If S* is close to 1 (epidemic at its start), then βS* > γ, and it's necessary to choose Δt very small to respect the CFL condition.

The explicit Euler scheme used is of order 1 in time, which means that the error decreases as O(Δt). So convergence is guaranteed if the scheme is stable and consistent with the initial differential equation system.

Implementing the simulation with the current method gives us solutions that we consider as reference solutions for this study. The following graph shows the evolution of the populations for the recorded disease over 100 days.

<figure><img src="image_4.png" /><figcaption>Figure 1: Evolution of populations (FDM)</figcaption></figure>
## 4. Markov chain approximation

Let E be the set of model states. Then we have:

$$E = \{S, I, R\}$$

In the rest of the study, we'll assume that a recovered individual now has permanent immunity to the disease, and can therefore no longer become susceptible.

### Formalization

* $S \to I$: individuals go from being “Susceptible”to “Infected”with a probability proportional to the number of infected individuals $\frac{\beta I}{N}$.

**NB:** The proportion is normalized here because we're looking for a probabilistic equivalence.

* $I \to R$: Individuals go from “Infected”to “Recovered”with a probability proportional to $\gamma$.

We consider the following graph:

$$1 - \frac{\beta I}{N} \subset S \xrightarrow{\frac{\beta I}{N}} I \xrightarrow{\gamma} R \supset 1$$

Let P be a transition matrix with $P_{ij}$ the transition probability from state i to j.

$$P = \begin{pmatrix} P(S \to S) & P(S \to I) & P(S \to R) \\ P(I \to S) & P(I \to I) & P(I \to R) \\ P(R \to S) & P(R \to I) & P(R \to R) \end{pmatrix}$$

where

$$P(S \to R) = P(I \to S) = P(R \to S) = P(R \to I) = 0, P(R \to R) = 1$$

So

$$P = \begin{pmatrix} P(S \to S) & P(S \to I) & 0 \\ 0 & P(I \to I) & P(I \to R) \\ 0 & 0 & 1 \end{pmatrix}$$

$$P = \begin{pmatrix} 1 - \frac{\beta I}{N} & \frac{\beta I}{N} & 0 \\ 0 & 1 - \gamma & \gamma \gamma 0 & 0 & 1 \end{pmatrix}$$

Let $\Delta I$ and $\Delta R$ be the respective numbers of new infections and new recoveries.

If $S_t$ is the number of Susceptibles at time t, then the number of newly infected would be the number of successes (k) in $S_t$ independent trials, where each trial has $\frac{\beta I}{N}$ as the probability of success.
Thus $\Delta I \sim \mathcal{B}(S_t, \beta \frac{I}{N})$ and $\Delta R \sim \mathcal{B}(I_t, \gamma)$

## Evolution of populations

Let $X_t = (S_t, I_t, R_t)^T$ be the “population”vector at time t. The time evolution is defined by:

$$X_{t+1} = P_t X_t$$

where $P_t$ is the transition matrix at time t.

So, knowing the populations at time t, we can easily get an estimate at time t+1.

Applying this method to the above problem, we obtain the following graph of temporal evolution

<figure><img src="image_5.png" /><figcaption>Figure 2: Evolution of populations (Markov)</figcaption></figure>

We note that the trend is generally captured by referring to the reference solution, apart from a few irregularities. This has no impact on the interpretation of the model and the decision to eradicate this disease. A detailed analysis of the results obtained is given below.

# 5. Neural Network approximation

Physics-Informed Neural Networks (PINNs) are an innovative approach to solving equations and systems of differential equations using neural networks.

## PINNs setup

The aim of PINNs is to approximate the functions $S(t)$, $I(t)$, and $R(t)$ by a neural network $u_\theta(t) = [S_\theta(t), I_\theta(t), R_\theta(t)]$, where $\theta$ represents the network parameters. The training aims to minimize the residuals of the differential equations and respect the initial conditions.
## Residuals of Differential Equations

The residuals for each equation are defined as follows

$$ \mathcal{R}_S(t, \theta) = \frac{\partial S_{\theta}}{\partial t} + \beta S_{\theta} I_{\theta}, $$

$$ \theta \mathcal{R}_I(t, \theta) = \frac{\partial I_{\theta}}{\partial t} - (\beta S_{\theta} I_{\theta} - \gamma I_{\theta}), $$

$$ \theta \mathcal{R}_R(t, \theta) = \frac{\partial R_{\theta}}{\partial t} - \gamma I_{\theta}. $$

The aim is to minimize these residuals to ensure that $u_{\theta}(t)$ satisfies the differential equations.

Initial conditions are imposed via a penalty in the loss function (2):

$$ \mathcal{L}_{\text{IC}} = \|S_{\theta}(0) - S_{0}\|^{2} + \|I_{\theta}(0) - I_{0}\|^{2} + \|R_{\theta}(0) - R_{0}\|^{2}. $$

## Total Loss Function

The total loss function combining residuals and initial conditions is given by :

$$ \mathcal{L}(\theta) = \mathcal{L}_{\text{ODE}} + \mathcal{L}_{\text{IC}}, $$

with

$$ \mathcal{L}_{\text{ODE}} = \int_{0}^{T} (\mathcal{R}_{S}^{2}(t, \theta) + \mathcal{R}_{I}^{2}(t, \theta) + \mathcal{R}_{R}^{2}(t, \theta)) dt. $$

In practice, this integral is approximated by a discrete sum over a set of collocation points (Ghosh) $\{t_i\}_{i=1}^N$:

$$ \mathcal{L}_{\text{ODE}} \approx \frac{1}{N} \sum_{i=1}^{N} (\mathcal{R}_{S}^{2}(t_i, \theta) + \mathcal{R}_{I}^{2}(t_i, \theta) + \mathcal{R}_{R}^{2}(t_i, \theta)). $$

Convergence is obtained if :

* $\mathcal{L}(\theta) \to 0$ during training,

* The network is flexible enough to represent the exact solution.

For this part, we use neural networks with 3 different numbers of iterations in order to assess performance and approximate the reference solution as closely as possible.

We also accompany these graphs of temporal evolutions with the loss curve to visualize how well the model performs and minimizes errors over the course of iterations. The network was implemented using the **PyTorch** (PyTorch), **numpy** and ** matplotlib** packages. The results obtained for different numbers of iterations are shown below:
<figure><img src="image_6.png" /><figcaption>Figure 3: PINNs with 1000 iterations</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Loss - Total training time=11.9s</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: PINNs with 5000 iterations</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Loss - Total training time=1min 47.6s</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: PINNs with 10000 iterations</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Loss - Total training time=2min 46.2s</figcaption></figure>
A quick glance at the graphs tells us that near-perfect convergence is reached after 10000 iterations, but an interpretation similar to that of the reference solution can also be derived from the solution with 5.000 iterations.

# 6. Discussion

The analysis of the simulation graphs provides valuable insight into the progression of disease X within the population over time. Upon closer inspection, it becomes evident that approximately 30 days into the simulation, there is a gradual reduction in the number of individuals infected with disease X, which continues steadily until the final 70 days under consideration. By this point, the population is predominantly composed of individuals who have recovered from the disease, with permanent immunity, alongside a minor fraction of individuals who remain susceptible and could potentially contract the disease. Notably, the highest recorded number of infected individuals peaks at around 800, which corresponds to approximately 32% of the total population. This observation strongly suggests that disease X spreads rapidly within the population, with an initial reproduction number $R_0 = \frac{\beta}{\gamma} = 2.99 > 1$, indicating that, on average, each infected individual is transmitting the disease to more than one other individual.

These simulation graphs exhibit typical epidemic dynamics under controlled conditions, where the number of susceptible individuals diminishes, the infected population rises sharply before gradually declining, and ultimately, the number of recovered individuals predominates in the long run. The early peak in infections followed by a swift decline in active cases underscores the rapid, yet controlled, spread of the disease within the simulated environment.

When using the finite-difference method as a benchmark, it becomes apparent that the proposed Markovian approach is highly effective in accurately capturing and modeling the progression of disease X in the specified region (denoted as Y). Despite some minor fluctuations in the data, these irregularities have no significant impact on the overall interpretation and analysis of the phenomenon. The Markovian model proves to be a reliable representation of the disease's dynamics over time.

Regarding the use of Physics-Informed Neural Networks (PINNs), it is noteworthy that even after as few as 5,000 iterations, the predictions generated by the neural network model closely follow the expected evolution of disease X. After 10,000 iterations, the PINN-based predictions align seamlessly with the results obtained from traditional numerical methods, such as the finite-difference method (MDF) and the Markovian model, demonstrating the ability of the neural network approach to effectively replicate the dynamics captured by conventional techniques.

In conclusion, both the Markovian approach and the PINNs-based methodology have proven to be effective and robust tools for modeling the spread of disease X in the population. Their ability to generate reliable predictions offers a solid foundation for the development of future preventive measures, helping to inform decision-making processes aimed at managing the spread of infectious diseases. The results of these models provide valuable insights that can be leveraged in real-world scenarios to mitigate the impact of the disease.
# 7. Conclusion

The study of the progression of disease X through various approaches, including the finite-difference method, the Markovian approach, and Physics-Informed Neural Networks (PINNs), has effectively demonstrated the utility and accuracy of these methods in modeling the phenomenon. The Markovian approach and PINNs stand out due to their ability to faithfully describe the evolution of the disease and provide predictions that align closely with the results obtained through traditional numerical methods. These tools prove to be highly relevant for understanding the dynamics of disease transmission and can be used as a basis for formulating appropriate preventive strategies. Their capacity to offer reliable insights into disease progression makes them invaluable for decision-making processes aimed at controlling and mitigating the impact of the disease.

# References

[1] Andersson, H. and Britton, T. (2000). Stochastic epidemic models and their statistical analysis. In *Lecture Notes in Statistics*, Springer Verlag, pages 11–18.

[2] Caterina Millevoi, D. P. and Ferronato, M. (2023). A physics-informed neural network approach for compartmental epidemiological models. In *axXiv math.NA*.

[Ghosh] Ghosh, K. Time series forecasting for epidemic prediction using pinn and lstm. In *Indian Institute of Science, Bengaluru*.

[Lapebie] Lapebie, E. Epidémiologie mathématique : 2. les modèles à compartiments. In https://apex-solutions.fr/epidemiologie-mathematique-2-les-modeles-a-compartiments/.

[PyTorch] PyTorch. Documentation. In https://pytorch.org/docs/stable/index.html.