Numerical Scheme for Compartmental Models: New Matlab Software Codes for Numerical Simulation

Introduction

With the help of a programming language that represents matrix and array mathematics directly, MATLAB combines a desktop environment tailored for iterative analysis and design processes.¹ The Live Editor for writing scripts that mix code, output, and formatted text in an executable notebook is part of it.¹ The Windows requirements are Windows 10 (version 20H2 or higher), Windows 11, Windows Server 2019, and Windows Server 2022. Many scientists and mathematicians choose to use this program since it can be accessed across all popular platforms, including Linux, Mac, and Windows, and because it can be used to explore, model, and analyze data.¹ Matlab software has been used to run numerical simulations of compartmental models in epidemiology.^2–11

There are several fundamental compartmental models described using differential equations. The basic ones include Susceptible - Infected (SI), Susceptible - Infected - Recovered (SIR),^12,13 Susceptible - Infected - Susceptible (SIS),¹⁴ Susceptible - Infected - Recovered - Vaccinated (SIRV),¹⁵ Susceptible - Exposed - Infected - Recovered (SEIR) models.⁴ The purpose of this study is to make public new Matlab codes that authors have been utilising in their work to aid researchers, especially students, who rely on deterministic or compartmental modeling of epidemiology in the numerical simulation of their research projects. This well-detailed code, in our opinion, might be extremely helpful to them as many of them struggle to do the numerical simulations due to the dearth of research that specifically tackles numerical simulation of deterministic models and also to provide users more freedom for coding in Matlab. Recently, researchers have started sharing their codes and providing detailed explanation on how to use them. To provide users extra coding freedom, Guo et al.¹⁶ presented newly developed visualization framework called OpenSeesPyView, which is a Python programming-based graphical user interface (GUI) for OpenSeesPy, a prevalent finite element solver in earthquake engineering. The R package ag5Tools was written by Brown et al.¹⁷ and offers a streamlined interface for downloading and retrieving AgERA5 data. With the help of the program, time-series data for groups of geographic points may be easily extracted and converted into a format that can be employed in statistical models used in agricultural research. The Rcall interface, developed by Egert and Kreutz,¹⁸ gives users access to a large range of techniques written in MATLAB and R. The program is MATLAB-based and offers direct access to R-based tools and methodologies, such as those found on Bioconductor or CRAN. ShinyGAStool, an open source tool created by Hoffmann et al.,¹⁹ allows users to easily execute a candidate gene association analysis from a web browser using huge datasets. The remaining section are group as follows: The method section, where we demonstrate how to use the Matlab software codes. we look at the implementation, operation and discussions and limitations. The last section is the concluding section. The software section has the alternative software codes in R.

Methods

Implementation

We demonstrate how to use the matlab software codes with the SEIR compartmental model depicted in Figure 1. The population is partitioned into four (4) compartments: Susceptible, exposed, infected and recovered. Individuals are recruited into the susceptible class at a rate $\mathrm{\Omega }$ and they die at a rate $\mu$. The transmission rate is $\beta$ and the recovery rate is $\gamma$. The rate at which exposed individuals become infectious is $\alpha$ and the disease-induced death rate is $\sigma$.

Figure 1. Model dynamics flowchart.

The model is described by the following ordinary differential equations.

(1)

$$\begin{array}{l}\frac{\mathit{dS}}{\mathit{dt}}=\mathrm{\Omega }-\beta \frac{\mathit{\text{IS}}}{N}-\mathit{\mu S},\\ \frac{\mathit{dE}}{\mathit{dt}}=\beta \frac{\mathit{\text{IS}}}{N}-\left(\alpha +\mu \right)E,\\ \frac{\mathit{dI}}{\mathit{dt}}=\mathit{\alpha E}-\left(\gamma +\sigma +\mu \right)I,\\ \frac{\mathit{dR}}{\mathit{dt}}=\mathit{\gamma I}-\mathit{\mu R}.\end{array}$$

with initial conditions $S\left(0\right)=S_0,E\left(0\right)=E_0,I\left(0\right)=I_0,R\left(0\right)=R_0$.

For the purposes of the simulations, the following parameter values are chosen and is given in Table 1. The initial conditions chosen are $S\left(0\right)=1000,E\left(0\right)=0,I\left(0\right)=1,R\left(0\right)=0$

Table 1. Parameter values and description.

Parameter	Description	Value	Source
$\mathrm{\Omega }$	Recruitment	29.08	2
$\beta$	Transmission rate	0.9	2
$\alpha$	infectiousness of the exposed individuals	0.3	assumed
$\mu$	natural rate of death	0.4252912 $\times {10}^{-4}$	2
$\gamma$	Recovery rate of infected individuals	0.1	assumed
$\sigma$	Disease-induced death rate	0.003286	assumed

Once you have your model and parameter values clearly defined, you can then open the matlab editor window which is shown in Figure 2. Input or copy the codes and paste at the new script.

Figure 2. Matlab Editor window.

The numerical Matlab software codes used for model (1) $\mathit{\text{function}}\left[t,S,E,I,R\right]=\mathit{\text{SEIRMODEL}}\left(\mathrm{\Omega },\beta ,\mu ,\sigma ,\alpha ,\gamma ,N,S_0,E_0,I_0,R_0,\mathit{\text{MaxTime}}\right)$

Operation

if nargin==0

$\mathrm{\Omega }=20;\beta =0.3;\mu =0.00004252912;\alpha =0.3;\sigma =0.003286;\gamma =0.1;$

$S_0=1000;E_0=0;I_0=1;R_0=0.0;$

MaxTime=120; % The MaxTime stands for the duration of the simulation

end

$S=S_0;E=E_0;I=I_0;R=R_0;$ N=S+E+I+R;

% The main iteration

options=odeset(’RelTol’,1e-5);

$\left[t,\mathit{pop}\right]=\mathit{ode}45\left(@\mathit{\text{Diff}}26,\left[0\mathit{\text{MaxTime}}\right],\left[\mathit{\text{SEIR}}\right],\mathit{\text{options}},\left[\mathrm{\Omega }\gamma \beta \mu \sigma \alpha N\right]\right);$

$S=\mathit{pop}\left(:,1\right);E=\mathit{pop}\left(:,2\right);I=\mathit{pop}\left(:,3\right);R=\mathit{pop}\left(:,4\right);$

% plots the graphs with scaled colours

figure(1)

Y=plot(t,S,’b.’);

legend(Y,’Susceptible’)

xlabel’Time (days)’

ylabel’S(t)’

figure(2)

T=plot(t,E,’.g’);

legend(T,’Exposed’)

xlabel’Time (days)’

ylabel’E(t)’

figure(3)

f=plot(t,I,’.m’);

legend(f,’Infected’)

xlabel’Time (days)’

ylabel’I(t)’

figure(4)

h=plot(t,R,’.g’);

legend(h,’Recovered’)

xlabel’Time (days)’

ylabel’R(t)’

% calculates the differential rates used in the integration.

$$\begin{array}{l}\mathit{\text{function}}\mathit{\text{dpop}}=\mathit{\text{Diff}}26\left(t,\mathit{pop},\mathit{\text{parameter}}\right)\\ \mathrm{\Omega }=\mathit{\text{parameter}}\left(1\right);\gamma =\mathit{\text{parameter}}\left(2\right);\beta =\mathit{\text{parameter}}\left(3\right);\mu =\mathit{\text{parameter}}\left(4\right);\\ \sigma =\mathit{\text{parameter}}\left(5\right);\alpha =\mathit{\text{parameter}}\left(6\right);N=\mathit{\text{parameter}}\left(7\right);\\ S=\mathit{pop}\left(1\right);E=\mathit{pop}\left(2\right);I=\mathit{pop}\left(3\right);R=\mathit{pop}\left(4\right);\\ \mathit{\text{dpop}}=\mathit{\text{zeros}}\left(4,1\right);\%\mathit{\text{These}}\mathit{are}\mathit{\text{the}}\mathit{\text{system}}{\left(1\right)}^{\prime }s\mathit{\text{equations}}\\ \mathit{\text{dpop}}\left(1\right)=\mathrm{\Omega }-\left({\beta }^{\ast }{I}^{\ast }S\right)./\left(N\right)-{\mu }^{\ast }S;\\ \mathit{\text{dpop}}\left(2\right)=\left({\beta }^{\ast }{I}^{\ast }S\right)./\left(N\right)-\left(\alpha +\mu \right)\ast E;\\ \mathit{\text{dpop}}\left(3\right)={\alpha }^{\ast }E-\left(\gamma +\mu +\sigma \right).\ast I;\\ \mathit{\text{dpop}}\left(4\right)={\gamma }^{\ast }I-{\mu }^{\ast }R;\end{array}$$

Upon running the codes, some of the simulation results are shown in Figures 3–6.

Figure 3. Exposed compartment.

Figure 4. Infected Compartment.

Figure 5. Recovered Compartment.

Figure 6. Susceptible compartment.

Discussion

In Figures 3–6, are reported, the numerical solutions of system (1) for a period of 120 days. These codes can be modified for any compartmental models. The ‘figure’ command produces the output given in Figures 3–6. The steps or procedures listed in the codes have to be followed carefully in order not to encounter errors. The parameters can be represented with letters for instance $\mathrm{\Omega }$ can be written in the codes as Omega as Matlab doesn’t recognize the parameters listed in the code. The code written in the editor window can be seen in Figure 7. The SEIR model is extended, and an alternative software, R codes has been provided at the appendix section. Using the same initial conditions and parameter values given in Table 1, the output figures for the R code are given by Figures 8–10. Users who cannot afford Matlab software can freely use the R software codes for the numerical simulation. The two software codes gives the same output results.

Figure 7. Matlab software codes in editor window.

Figure 8. Susceptible compartment using R.

Figure 9. Exposed compartment using R.

Figure 10. Infected compartment using R.

Conclusion

This work seeks to introduce new matlab software codes purposely for numerical simulations of classical compartmental models which can run on any version of Matlab. The intended targets are researchers and students who uses Matlab for their analysis. This codes can be edited/modify to suit any deterministic models and any desire output required. The SEIR deterministic model was used to give a much insight about the codes. Alternatively, a deterministic SEIR codes written in R software is provided for those who wants to use freely available software. Despite the limitations of the R software, the deterministic model implemented in the R code can still be a useful tool for understanding the basic dynamics of disease transmission.