# Application of neural-network hybrid models in estimating the infection
functions of nonlinear epidemic models

Chentong Li  Guangdong Key Laboratory of Modern Control Technology, Institute
of Intelligent Manufacturing, Guangdong Academy of Science, Guangzhou,
Guangdong 510070, P.R.China  Changsheng Zhou  School of Mathematics and
Information science, Guangzhou University, Guangzhou, Guangdong 510006,
P.R.China  Junmin Liu  School of Mathematics and Statistics, Xi’an Jiaotong
University, Xi’an, Shaanxi 710049, P.R.China  Yao Rong  College of Engineering
Physics, Shenzhen Technology University, Shenzhen, Guangdong 518118, P.R.China

###### Abstract

Hybrid neural-network models combine the advantages of a neural network’s
fitting functionality with differential equation models to reflect actual
physical processes and are widely used in analyzing time-series data. Most
related studies have focused on linear hybrid models, but only a few have
examined nonlinear problems. In this work, we use a hybrid nonlinear epidemic
neural network as the entry point to study its power in predicting the correct
infection function of an epidemic model. To achieve this goal, we combine the
bifurcation theory of the nonlinear differential model with the mean-squared
error loss and design a novel loss function to ensure model trainability.
Furthermore, we find the unique existence conditions supporting ordinary
differential equations to estimate the correct infection function. Using the
Runge–Kutta method, we perform numerical experiments on our proposed model and
verify its soundness. We also apply it to real COVID-19 data to accurately
discover the change law of its infectivity.

Key Words: Differential equations; Epidemic model; Hybrid model; Neural
network.

##  1 Introduction

Forecasting an epidemic’s infection number and rate is a key objective in the
study of infectious diseases. The ordinary differential equation model of
infectious disease is based on known disease properties and has long been used
to analyze the dynamic behaviors of infection and to provide forecasts [4, 2,
9]. All infectious disease models can be divided into two parts: infection
terminology and susceptibility [2, 13]. The infected terms are used to
illustrate the change laws of the number of infected persons at different
periods (e.g., incubation, infection, and isolation). Most epidemic models use
nonlinear infection functions. However, nonlinearity quickly spawns from the
many combinations of infection and susceptibility terms [13]. To deal with the
difficulties of nonlinearity, researchers usually use the bifurcation theory
of nonlinear equations to make estimations; then, they improve upon them with
expert knowledge [13].

To implement nonlinear epidemic models, the basic reproduction bifurcation
number, $R_{0}$, comprises the parameters and functions of the model [23].
This number is used to describe the expected number of infected cases based on
one original infected person [24]. Used as a biological descriptor, $R_{0}$ is
the threshold value used to determine whether the disease may vanish. In many
epidemic models, when $R_{0}$ is larger than $1$, the epidemic will persist
forever; if it is smaller than $1$, the number of infected individuals will
approach $0$ over time [24]. This persistence description is called “forward
bifurcation” [14] and is the key output of epidemic modeling.

Although nonlinear epidemic model analysis has come a long way, selecting the
correct infection function (also known as the infection incidence) remains
prohibitively difficult [29, 20, 1, 3, 19, 26]. The functions not only pertain
to the number of susceptible and infected individuals, but they also involve
environmental and human behavior factors. For example, the periodical
temperature changes in temperate zones may influence the activity of viral
proteins, causing infectability to reflect a periodical property [29, 20].
Additionally, when the infected number grows, the growth rate of infectability
may become smaller, depending on human behaviors. Mathematically, the Holling
type function is applied to the infection function to account for such
behavior. However, the many combinations of factors create considerable
uncertainty [19, 26]. Thus, any inaccuracy in the choice of infection function
directly influences the precision of model prediction .

Some statistical methods have been applied to real epidemic data using
Akaike’s information criterion (AIC) [27, 28] or its second-order estimate
(AICC) [10]. Using these, one can ostensibly find the most suitable function
from a group of candidates. In many cases, however, relying on the available
data is insufficient, and expert knowledge of the environment, human behavior,
and the properties of the viral strain is also needed. Therefore, better
methods are necessary.

Notably, improvements to deep learning [11] and related neural-network models
have greatly enhanced the predictability of nonlinear functions throughout
modern industry. A neural network is modeled after the neuromorphology of the
human brain using network to construct input layers, output layers, neurons,
and activation functions (synapses) at each layer. As with many physical
functions, it sacrifices explainability for precision. Hence, at its heart is
a black-box function that is trained to “think” in different dimensions [7].
According to the universal approximation theory [8], a neural network can
approximate any functions given the ssoftwareuitable activation functions.
Therefore, if we can build a neural network to implement the differential
equations of an infectious disease model, infection model prediction should
become much more accurate.

As introduced by Psichogios and Ungar [16], combining capabilities in this way
results in a “hybrid” neural-network model. The authors built a hybrid model
to forecast biological phenomena. Later, [5, 12] used them to predict the
change laws of time-series data. Recently, the authors of [18, 17, 25] had
breakthrough success in forecasting biological processes. Notably, biological
systems contain both first-principle information (i.e., biophysical processes)
and unknown information, which accounts for the complexity of a biophysical
system. Hence, the problem now reduces to properly fitting the data to the
biophysical process.

Most hybrid neural-network studies have leveraged the linear-first principle
to effectively side-step overly complex nonlinear preparation methods. Thus,
questions have been raised about whether the hybrid functions can be
adequately trained with limited data and whether the estimates are actually
reliable. In this work, based on the standard theorem of forwarding
bifurcation, a loss function is created that overcomes the untrainability of
the nonlinear epidemic model. Furthermore, based on the unique conditions of
ordinary differential equations and the basic properties of epidemic models,
we provide the conditions needed to predict the correct infection function.

The remainder of this work is organized as follows. In the second section, we
provide complete descriptions of the nonlinear epidemic system and introduce
the basic conditions of ordinary differential epidemic modeling. In the third
section, the two theorems that can solve the trainable and estimable problems
of the neural network are explained. In section four, based on the standard
Runge–Kutta method of numerically solving ordinary differential equations, we
introduce a method of properly training the hybrid neural network. In the
fifth section, the results of four numerical experiments are reported to
demonstrate the efficacy of the proposed model in estimating periodical and
Holling-type infection functions. Then, in the sixth section, we report on the
application of the proposed method to real data. In the last section, we
conclude this work and discuss future opportunities.

##  2 Hybrid epidemic model properties and conditions

The hybrid nonlinear epidemic neural network model is written as follows:

$x^{\prime}=g(x,f),$ (1)

where $x=(x_{S},x_{I})$ and $x_{S}\in R^{m_{1}}$ are the variables of a
susceptible person, and $x_{I}\in R^{m_{2}}$ are those of an infected person.
Function $f$ maps time $t$ or variables $x$ into $R$. Function
$g=(g_{S},g_{I})$, where $g_{S}$ is the function mapping of $x_{S}$ and the
output of function $f$ from $R^{m_{1}+1}$ to $R^{m_{1}}$. Additionally,
$g_{I}$ is the function mapping of $x_{I}$ and the output of function $f$ from
$R^{m_{2}+1}$ to $R^{m_{2}}$. Furthermore, we assume that function $g$
satisfies the following three conditions:

(C1) $\forall x_{1}\neq x_{2}$, there exist a positive constant, $L_{1}$, such
that $\|g(x_{1},y)-g(x_{2},y)\|<L_{1}\|x_{1}-x_{2}\|$;

(C2) $\forall y_{1}\neq y_{2}$, there exist a positive constant, $L_{2}$, such
that $\|g(x,y_{1})-g(x,y_{2})\|<L_{2}\|y_{1}-y_{2}\|$;

(C3) if $x_{I}\neq 0$, then $\forall y_{1}\neq y_{2}$, and the in-equation
$\|g_{I}(x,y_{1})-g_{I}(x,y_{2})\|>0$ hold.

Condition (C1) is the famous Lipschitz condition of the ordinary differential
equation and is sufficient for the equation’s unique existence [22]. Condition
(C2) is the Lipschitz condition for the function, which, when used as the
parameter, accommodates the unique existence of its corresponding ordinary
differential equation. In the proofs shown in the next section, we replace
$L_{1}$ and $L_{2}$ in (C1) and (C2), respectively, with
$L=\max(L_{1},L_{2})$. Condition (C3) guarantees that, when $x_{I}\neq 0$, map
$g$ is the injective map for $y$. These three conditions are used in the
upcoming proofs.

Let $x(x_{0},f,t)=(x_{S}(x_{0},f,t),x_{I}(x_{0},f,t))$ be the solution to
Equation (1) with respect to function $f$ and the initial condition, $x_{0}$.
Then, based on the standard of nonlinear epidemic models [13], we can assume
that two more conditions are satisfied by the system of Equation (1).

Condition (C4), without $x_{I}$ in system $x_{S}^{\prime}=g_{S}(x_{S},f)$, has
the unique positive equilibrium, $S^{*}$.

Condition (C5) allows vector $0\in R^{m_{2}}$ to become the zero vector, and
when $x_{I}=0$, $g_{I}(x,f)=0$.

Conditions (C4) and (C5) are built upon the properties of infectious diseases.
If no disease exists, the system only contains susceptible individuals, and
the system has a unique positive state. This, the without-disease system will
remain in a stable no-disease state until a new infected person enters the
system. Therefore, the system of Equation (1) has a unique disease-free
equilibrium, $x_{DFE}=(S^{*},0)$.

(a) Forward bifurcation diagram. (b) Solutions under different conditions.

Figure 1: Forward bifurcation and solutions under different bifurcation
conditions.

Based on the standard forward bifurcation theory of the nonlinear epidemic
model [14], there exists a threshold value, $R_{0}>0$, the basic reproductive
number, which determines whether a disease outbreak occurs. This value applies
to most epidemic model s. Therefore, we assume that the system of Equation (1)
satisfies the following forward bifurcation condition:

In the (C6) condition, there exists a function, $R_{0}^{f}$, that maps the
parameters and function $f$ of the system of Equation (1) to a positive real
number, $R_{0}$. This map is continuous, and when $R_{0}<1$, the disease-free
equilibrium, $x_{DFE}$, of the ordinary differential equation system, Equation
(1), is globally stable; when $R_{0}>1$, the disease-free equilibrium is
broken.

Condition (C6) represents the standard result of most epidemic systems. Hence,
research on nonlinear epidemic models seek to find the basic reproductive
number [29, 23, 24]. Figure 1 illustrates a simple forward bifurcation
epidemic system. Figure 1(a) shows that when $R_{0}$ is smaller than $1$, the
disease-free equilibrium is stable, and when this number is larger than $1$,
the disease-free equilibrium becomes unstable. Meanwhile, Figure 1(b) shows
the solution to $x_{I}$ under different bifurcation conditions. When
$R_{0}>1$, the solution appears as the black line in Figure 1(b), and when
$R_{0}<1$, the solution appears as the blue line in Figure 1(b). The above six
conditions are the basic properties satisfied by most epidemic models.

##  3 Hybrid model trainability and estimability

To finish the proofs, we first introduce some more notations. Let the neural-
network model,
$f_{\theta}(x)=\phi_{n}\circ\cdots\circ\phi_{2}\circ\phi_{1}(x)$, where
$\phi_{i}(x)=\sigma(w_{i}x+b_{i})$ is the linear combination of the matrix,
$w_{i}$, and the bias vector, $b_{i}$, with the activation function, $\sigma$.
The constant, $n>0$, is the number of layers in the neural network,
$f_{\theta}$, and $\theta$ is the parameter set, including the weights in
$w_{i}$ and $b_{i}$. All activation functions used in this study consist of
elemental functions. We also note the norm,
$\|(x_{1}(t),x_{2}(t),\cdots,x_{m}(t))\|=\sum_{i=1}^{m}\int_{0}^{T}x_{i}(t)^{2}dt$,
where $(0,T)$ is the collected area of the data. Then, the following Lemma can
be directly obtained:

######  Lemma 1.

The neural-network function, $f_{\theta}:R^{n}\to R^{m}$, is continuous in
$x\in(0,T)^{n}$ with the bounded weight parameter, $\theta=(w,b)$, and
elemental activation functions $\sigma$.

This can be proven directly by using the continuities of the elemental
functions.

When estimating function $f$, the most common data type is the number of
infected individuals. Thus, we use the distance between the terms of infected
individuals of Equation (1) and the data of the infected individuals as the
loss function. The trivial formula of the loss can be written as
$L_{\theta}=\|x_{I}(x_{0},f_{\theta},t)-D(t)\|$, where $D(t)$ are the data.
However, this formula is not suitable for Equation (1). Based on the following
theorem, we prove that $f_{\theta}$ cannot be trained from data $D(t)$.

######  Theorem 1 (Vanishing gradient).

The system of Equation (1) satisfies conditions (C1)–(C6) if there exists a
parameter set, $A\neq\emptyset$, which is a subset of $R^{n}$, and
$\forall\theta\in A$, the basic reproductive number satisfies
$R^{f}_{0}(f_{\theta})<1$. Then, $\forall\epsilon>0$ and $D(t)$, there exists
a positive constant, $\delta$, positive initial condition $x_{0}$, and
parameters $\theta$, such that the gradient of the loss function,
$L_{\theta}=\|x_{I}(x_{0},f_{\theta},t)-D(t)\|$, with respect to the neural-
network parameters, satisfies $\nabla_{\theta}L_{\theta}<\epsilon$.

###### Proof.

When $\theta\in A$, by (C6), the disease-free equilibrium,
$x_{DFE}=(S^{*},0)$, is stable. Thus it follows that $\forall\epsilon>0$,
there exists a constant, $\delta>0$, for all initial conditions,
$\|x_{0}-x_{DFE}\|<\delta$, such that, for every $t\geq 0$,
$\|x_{I}(x_{0},f_{\theta},t)-0\|=\|x_{I}(x_{0},f_{\theta},t)\|<\epsilon/2$.
Meanwhile, for another parameter $\hat{\theta}\neq\theta$ and $\hat{\theta}\in
A$, we could find a $\hat{\delta}>0$, such that for all initial conditions
$\|x_{0}-x_{DFE}\|<\hat{\delta}$, we have
$\|x_{I}(x_{0},f_{\hat{\theta}},t)\|<\epsilon/2$. Then, we let
$\delta^{\prime}=\min(\delta,\hat{\delta})$, and
$\|x_{0}-x_{DFE}\|<\delta^{\prime}$. Thus, we obtain the following equations:

$\displaystyle L_{\theta}-L_{\hat{\theta}}$
$\displaystyle\leq\|x_{I}(x_{0},f_{\theta},t)-D(t)\|-\|x_{I}(x_{0},f_{\hat{\theta}},t)-D(t)\|$
,
$\displaystyle\leq(\|D(t)\|+\|x_{I}(x_{0},f_{\theta},t)\|)-(\|D(t)\|-\|x_{I}(x_{0},f_{\hat{\theta}},t)\|)$
, $\displaystyle\leq(\|D(t)\|+\epsilon/2)-(\|D(t)\|-\epsilon/2)$
$\displaystyle=\epsilon.$

By the continuity of $f_{\theta}$ (Lemma 1) and $R_{0}(x)$, there exists a
positive constant, $a$, such that the small ball, $B(\theta,a)\in A$.
Therefore, by the above equations, $\forall\hat{\theta}\in B(\theta,a)$, we
have $|L_{\theta}-L_{\hat{\theta}}|<\epsilon$, which follows
$\nabla_{\theta}L_{\theta}<\epsilon$. ∎

The blue line of Figure 1(b) shows one special solution of the epidemic model
when $R_{0}<1$, in which the solution approaches zero very fast. The same is
true for the other solution of the same model, making the difference between
them negligible. This causes the gradient of the loss function to become
rather small, which is intuitive.

In the real world, almost all detected diseases have basic reproductive
numbers that are larger than one. However, in the training process, the
trained function, $f_{\theta}$, may appear in an area that makes the number
smaller than one. Therefore, based on the above theorem, the hybrid model is
untrainable. For this reason, we introduce a loss function that combines the
standard loss with the basic reproductive number, $R_{0}$, which is the
bifurcation parameter of the nonlinear epidemic system, Equation (1). This
makes the loss function non-zero when $R_{0}$ is smaller than $1$.

$L_{\theta}=\|x_{I}(x_{0},f_{\theta},t)-D(t)\|+\alpha\max(1-R^{f}_{0}(f_{\theta}),0),$
(2)

where $\alpha>0$ is the hyper-parameter determined by prior knowledge. In this
loss function, when $R_{0}<1$, according to Theorem 1, loss
$\|x_{I}(x_{0},f_{\theta},t)-D(t)\|$ may be very small. However, when the
value of $\max(1-R^{f}_{0}(f_{\theta}),0)$ is larger than zero, and its
minimal value is zero, the gradient of that part, hence the entire loss
function, becomes non-zero. When $R_{0}>1$, the loss function becomes standard
again. For other hybrid nonlinear models with different bifurcation
parameters, we can construct suitable bifurcation items in the loss function
to ensure trainability.

From the above analysis, we have obtained a suitable loss function that allows
the hybrid nonlinear neural-network model to be trained. However, the question
remains of whether the infection function can be learned from real data.
Therefore, we introduce the following loss function:

$L_{\theta}=\|x_{I}(x_{0},f_{\theta},t)-x_{I}(x_{0},f^{*},t)\|+\alpha\max(1-R^{f}_{0}(f_{\theta}),0).$
(3)

We must find whether the real-world infection function, $f^{*}$, can be fitted
by neural-network $f_{\theta}$. Using the following theorem, we prove that the
system of Equation (1) fits function $f^{*}$ using the loss function of
Equation (3).

######  Theorem 2 (Existence).

The system of Equation (1) satisfies Conditions (C1)–(C6) if
$R^{f}_{0}(f^{*})>1$. Thus, $\forall\epsilon>0$, $f_{\theta}$ satisfies
$\|f_{\theta}-f^{*}\|<\epsilon$ if and only if there exists a parameter set,
$\theta$, such that the loss function, $L_{\theta}$, in Equation (3) satisfies
$L_{\theta}<\epsilon$.

###### Proof.

The solution to the system of Equation (1) with an initial time of zero and
$x_{0}$ as the initial condition is written as $x(x_{0},f_{\theta},\Delta
t)=x_{0}+\int_{0}^{\Delta t}g(x,f_{\theta})dt$. Therefore, we let
$x^{\theta}_{1}=x_{0}+\int_{0}^{\Delta t}g(x,f_{\theta})dt$,
$x^{*}_{1}=x_{0}+\int_{0}^{\Delta t}g(x,f^{*})dt$, and
$x^{\theta}_{n+1}=x^{\theta}_{n}+\int_{n\Delta t}^{(n+1)\Delta
t}g(x,f_{\theta})dt$, $x^{*}_{n+1}=x^{*}_{n}+\int_{n\Delta t}^{(n+1)\Delta
t}g(x,f^{*})dt$ for $n=1,2,..,N$, where $N=\lfloor T/\Delta t\rfloor$.

By the definition of $\|\cdot\|$, function
$\|x(x_{0},f_{\theta},t)-x(x_{0},f^{*},t)\|$ is monotone increasing with
respect to $t$. Then, by Conditions (C1) and (C2), the following system of
equations holds:

$\displaystyle\|x^{\theta}_{1}-x^{*}_{1}\|$ $\displaystyle\leq\int_{0}^{\Delta
t}\|g(x,f_{\theta})-g(x,f^{*})\|dt$ $\displaystyle\leq\Delta
tL(\|f_{\theta}-f^{*}\|+\|x^{\theta}_{1}-x^{*}_{1}\|)$
$\displaystyle\leq\Delta tL(\epsilon+\|x^{\theta}_{1}-x^{*}_{1}\|).$

Let $\Delta t=\frac{1}{2L}$ be the time interval, so that we have
$\|x^{\theta}_{1}-x^{*}_{1}\|\leq\frac{\Delta tL}{1-\Delta
tL}\epsilon=\epsilon$. Similarly, the following hold:
$\|x^{\theta}_{2}-x^{*}_{2}\|\leq\frac{\Delta tL+1}{1-\Delta
tL}\epsilon=3\epsilon$, and $\|x^{\theta}_{n}-x^{*}_{n}\|\leq\frac{\Delta
tL+n-1}{1-\Delta tL}\epsilon=(2n-1)\epsilon$. Thus,

$\displaystyle L_{\theta}$
$\displaystyle=\|x_{I}(x_{0},f_{\theta},t)–x_{I}(x_{0},f^{*},t)\|\leq\sum_{n=1}^{N}\|x^{\theta}_{n}-x^{*}_{n}\|$
$\displaystyle\leq N^{2}\epsilon\leq(2TL+1)^{2}\epsilon.$

Let $\epsilon’=(2TL+1)^{-2}\epsilon$, so that we get
$L_{\theta}<\epsilon^{\prime}$.

On the other hand, supposing there exists a parameter set, $\theta^{\prime}$,
and a positive constant, $\delta$, such that
$\|f_{\theta^{\prime}}-f^{*}\|>\delta$ and $L_{\theta^{\prime}}<\epsilon$, it
follows that if $\|f_{\theta}-f^{*}\|<\epsilon$, then
$L_{\theta}<\epsilon^{\prime}$ by the first part of that proof, and
$\|f_{\theta^{\prime}}-f_{\theta}\|>\delta-\epsilon>0$. Thus, there exists a
time, $t_{c}\in(0,T)$, $\forall t\in[0,t_{c})$,
$f_{\theta^{\prime}}(t)=f_{\theta}(t)$ and $f_{\theta^{\prime}}(t_{c})\neq
f_{\theta}(t_{c})$. Therefore,

$\displaystyle\|x_{I}(x_{0},f_{\theta},t)-x_{I}(x_{0},f_{\theta^{\prime}},t)\|$
$\displaystyle\geq\int_{0}^{t_{c}+\delta
t}(x_{I}(x_{0},f_{\theta},t)-x_{I}(x_{0},f_{\theta^{\prime}},t))^{2}dt$
$\displaystyle=\int_{t_{c}}^{t_{c}+\delta
t}(\int_{t_{c}}^{t}g_{I}(x,f_{\theta})-g_{I}(x,f_{\theta^{\prime}})ds)^{2}dt.$

According to $\|g_{I}(x,f_{\theta})-g_{I}(x,f_{\theta^{\prime}})\|>0$
(Condition (C3)) and the continuity of function $g(x,y)$, there exists a
small-enough $\delta t$ that, $\forall t\in(t,t+\delta t)$,
$\|g_{I}(x(t),f_{\theta})-g_{I}(x(t),f_{\theta^{\prime}})\|>0$, from which it
follows that constant $\delta^{\prime}>0$ exists, such that
$\|g_{I}(x(t),f_{\theta})-g_{I}(x(t),f_{\theta^{\prime}})\|>\delta^{\prime}$.
Hence, we have
$\|x_{I}(x_{0},f_{\theta},t)-x_{I}(x_{0},f_{\theta^{\prime}},t)\|>\delta^{\prime
2}\delta t^{3}/3:=\hat{\delta}>0$. Thus,
$L_{\theta^{\prime}}>\hat{\delta}-\epsilon^{\prime}>0$, which leads to a
contradiction. ∎

From the above theorem, if the epidemic hybrid model satisfies the conditions
listed in the last section , the form of the infection function in the hybrid
model can be estimated using real data. In the next two sections, we explain
how to use a suitable numerical method to estimate the infection function.

##  4 Numerical methods

Figure 2: Calculation process.

Euler and Runge–Kutta methods are the most effective and most widely used for
solving ordinary differential equations [6]. The Runge–Kutta method derives
from the integral form of the ordinary differential equation and is used to
calculate the numerical solution at each time point. The Euler method derives
from the discrete form of differentiation, and it is less accurate. However,
it requires fewer computing resources.

In this work, given a known infection function with small perturbations, the
Euler method is used to generate training data. During training, the fourth-
order Runge–Kutta method with fixed steps is used to calculate the loss
between generated data and the hybrid model. Figure 2 shows the calculation
diagram of the numerical method. At each time step of the numerical
calculation, the neural-network input is $f_{in}(i\Delta t)$,
$i=1,2,\cdots,M$, where the $f_{in}(i\Delta t)=i\Delta t$, or $f_{in}(i\Delta
t)=x(i\Delta t)$, and the output is the value of the estimated function at
that time. Then, the numerical solution of the hybrid model at each step is
used to fit the generated data. In the loss function of Equation (2), the norm
of $\|\cdot\|$ is difficult to calculate directly. Thus , we use the discrete
form of $\|\cdot\|$ (mean-squared error (MSE) loss) to calculate its norm
[21]. By combining the MSE loss and the bifurcation items, we obtain the loss
function used for numerical calculation:

$L_{\theta}=\sum_{i=1}^{M}(y_{i,\theta}-D_{i})^{2}/M+\alpha\max(1-R^{f}_{0}(f_{\theta}),0),$
(4)

where $M=\lfloor T/dt\rfloor$, $dt$ is the length of the time step, and $T$ is
the length of the sample collected area.
$y_{i,\theta}=x_{I}(x_{0},f_{\theta},idt),andi=1,2,\cdots,M$ is the numerical
solution of the infected term of the system of Equation (1) at time $idt$.
$D_{i}$ is the collected data at time $idt$. The initial condition applies
$x_{0}$ number of susceptible and infected individuals. The number of infected
individuals reflects the initial number of cases, and the number of
susceptible persons reflects the entire population minus infected cases.

##  5 Numerical experiments

In this section, we report the results of four numerical experiments to
demonstrate the efficacy of the proposed hybrid neural-network model when
applied to estimating periodic and Holling type infection functions. All
experiments were performed on an Intel Core i7 -10700 @2.90-GHz computer with
16-GB RAM. All codes were written for the Windows x64 operating system using
Pytorch ($>=1.4.0$) and Numpy ($>=1.18.5$) Python packages. Source code can be
found at https://github.com/ChentongLi/Inf_Estimation_pytorch.

(a) Hybrid model fitted result. (b) Neural network fitted result. (c)
$|\beta(t)-\beta_{\theta}(t)|$. (d) Loss with respect to iterations.

Figure 3: Fitted results and susceptible–infected–susceptible training process
with a simply periodic function: (a) Fitted result of the hybrid model of
Equation (5) using generated data; (b) Fitted result of the neural network,
$\beta_{\theta}$, of the model of Equation (5) using the real function,
$\beta_{1}$, in the generation model; (c) Absolute value of errors between
$\beta_{\theta}$ and $\beta$; and (d) Change of loss function value over
iterations.

(a) Hybrid model fitted result. (b) Neural network fitted result. (c)
$|\beta(t)-\beta_{\theta}(t)|$. (d) Loss with respect to iterations.

Figure 4: Fitted results and training process of the
susceptible–infected–susceptible model with complex periodic function: (a)
Fitted result of the hybrid model of Equation (5) with the generated data; (b)
Fitted result of neural-network $\beta_{\theta}$ of the model of Equation (5)
with function $\beta_{2}$ in the generation model; (c) Absolute value of
errors between $\beta_{\theta}$ and $\beta$; and (d) Change of loss function
value over iterations.

The first and second numerical experiments were based on the standard
susceptible –infected–susceptible (SIS) model with birth and death data [13]:

$\left\\{\begin{aligned} \frac{dS}{dt}&=\lambda-\beta(t)SI+\gamma I-dS,\\\
\frac{dI}{dt}&=\beta(t)SI-\gamma I-dI,\end{aligned}\right.$ (5)

where the $S$ is the number of susceptible individuals, and $I$ is the number
of infected individuals at time $t$. Parameter $\lambda$ is the constant birth
rate, $\gamma$ is the recovery rate, and $d$ is the death rate. Function
$\beta(t)$ is the periodic function of time $t$ and is the infection function
about the infection incidence. In the first experiment, we used the simply
periodic function, $\beta_{1}(t)=0.5\Big{(}\cos(t/11)+1\Big{)}$, as the
infection function to generate data $D(t)$, and we used the hybrid model with
the neural network as the $\beta(t)$ to fit these data. In the second
experiment, the complex periodic function,

$\beta_{2}(t)=0.125\Big{(}-\cos(t/11)-\sin(t/11)-2.5\cos(2t/11)+0.5\sin(2t/11)\Big{)}+0.5,$

was used as the infection function. In the hybrid model of the second
experiment, the neural network was used to fit the infection function.
Equation (5) is a special case of the system of Equation (1), with $S=x_{S}$
and $I=x_{I}$. In these two experiments, the neural-network models were the
same. A three-layer neural network was used to equip the hybrid model with
fully connected layers of $1\times 16$, $16\times 16$, and $16\times 1$
matrices and their related biases. Both activation functions in the first and
second layers were the same: $\sin(ax)^{2}/a$ with $a=0.1$ [30]. The $R_{0}$
function used in the loss function of Equation (2) was numerically written as
$\frac{\sum_{i=1}^{M}\beta_{\theta}(i\Delta t)}{M(\gamma+d)}$ [23], where
$\beta_{\theta}(t)$ is the fitted function of the neural network, and the
constant, $M$, is the length of the time step. For hyperparameter $\alpha$ in
Equation (2), we chose $300(\gamma+d)$.

Figures 3 and 4 illustrate the results of the first experiment regarding
simple and complex periodic functions. During the training process of the
first experiment, the initial learning rate was set to $1\times 10^{-2}$ and
changed to $3.5\times 10^{-3}$ and $1.225\times 10^{-3}$ at the $250$th and
$1,500$th iterations, respectively. During the training process of the second
experiment, the initial learning rate was set to $1\times 10^{-2}$ and changed
to $3\times 10^{-3}$ and $9\times 10^{-4}$ at the $200$th and $1,500$th
iterations, respectively. In figures 3 and 4, subfigures (a), (b), (c), and
(d) present the respective fitted results of the hybrid model, the neural
network, their errors, and the training process. From these two experiments,
based on the loss function of Equation (4), we found that for the samples
collected (blue points), the hybrid model fit the data well, and the neural
network fit the infection function very well, which verifies the results of
Theorem 2. However, in the area lacking samples, the neural network did not
fit the infection function well. Comparing the results of Figure 3 with Figure
4, based on the same neural network model with the activation function,
$\sin(ax)^{2}/a$, we can see that the complex periodic method produced better
accuracy. During training, the simple periodic method had more vibrations.

(a) Hybrid model fitted result. (b) Neural network fitted result. (c)
$|\beta(S)-\beta_{\theta}(S)|$. (d) Loss with respect to iterations.

Figure 5: Fitted results and training process of the
susceptible–infected–recovered model with infection function $\beta(S)$: (a)
Fitted result of the hybrid model of Equation (6) with the generated data; (b)
Fitted result of neural-network $\beta_{\theta}$ of model of Equation (6) with
function $\beta(S)$ in the generation model; (c) Absolute value of errors
between $\beta_{\theta}$ and $\beta$; and (d) Change of loss function value
over iterations.

(a) Hybrid model fitted result. (b) Neural network fitted result. (c)
$|\beta(I)-\beta_{\theta}(I)|$. (d) Loss with respect to iterations.

Figure 6: Fitted results and training process of the
susceptible–infected–recovered model with infection function $\beta(I)$: (a)
Fitted result of the hybrid model of Equation (7) with the generated data; (b)
Fitted result of neural-network $\beta_{\theta}$ of the model of Equation (7)
with function $\beta(I)$ in the generation model; (c) Absolute value of errors
between $\beta_{\theta}$ and $\beta$; and (d) Change of loss function value
over iterations.

The third and fourth experiments were based on the susceptible–infected–
recovered (SIR) model with different infection functions correlated with
variables $S$ and $I$, which can be written as $\beta(S)$ and $\beta(I)$ in
this case. The model about the $\beta(S)$ is listed as

$\left\\{\begin{aligned} \frac{dS}{dt}&=-\beta(S)I,\\\
\frac{dI}{dt}&=\beta(S)I-\gamma I,\\\ \frac{dR}{dt}&=\gamma
I,\end{aligned}\right.$ (6)

where $S$ and $I$ are same as those in the SIS model, and variable $R$ is the
number of recovered individuals at time $t$. Parameter $\gamma$ is the
recovery rate. This model is used to generate data and construct the hybrid
model for the third experiment. Similarly, the model about $\beta(I)$ is
listed as

$\left\\{\begin{aligned} \frac{dS}{dt}&=-\beta(I)S,\\\
\frac{dI}{dt}&=\beta(I)S-\gamma I,\\\ \frac{dR}{dt}&=\gamma
I.\end{aligned}\right.$ (7)

This model is used in the fourth experiment. Infection incidences
$\beta(S)=0.2S/(S+20)$ and $\beta(I)=0.2I/(I+20)$ appear in the third and
fourth experiments, respectively. These two SIR models are equivalent to the
previous ones and represent special cases of the system of Equation (1) with
$I=x_{I}$ and $S=x_{S}$. In both experiments, we used a neural-network model
with three fully connected layers to fit the infection functions, and the
activation functions in the first and second layers were the same: $\tanh(x)$.
The three layers were constructed with $1\times 16$, $16\times 16$, and
$16\times 1$ matrices and their related biases. In the loss function of
Equation (4), the $R_{0}$ for the $\beta(S)$ system was
$\frac{\beta_{\theta}(S_{0})}{\gamma}$, and for the $\beta(I)$ system, it was
$\frac{10(\beta_{\theta}(0.1)-\beta_{\theta}(0))S_{0}}{\gamma}$, where $S_{0}$
was the initial condition of the hybrid model, and function $\beta_{\theta}$
was the neural-network model. The hyperparameter, $\alpha$, was set to
$300\gamma$.

Figures 5 and 6 respectively illustrate the results of the third experiment
about the infection function, $\beta(S)$, and the fourth experiment about the
infection function, $\beta(I)$. During the training process of the third
experiment, the initial learning rate was set to $1\times 10^{-2}$ and changed
to $4\times 10^{-3}$ and $1.6\times 10^{-3}$ at the $250$th and $3,000$th
iterations, respectively. During the training process of the fourth
experiment, the initial learning rate was set to $4\times 10^{-2}$ and changed
to $1\times 10^{-2}$, $2.5\times 10^{-3}$, and $6.25\times 10^{-4}$ at the
$250$th, $1,500$th, and $3,000$th iterations, respectively. In Figures 5 and
6, subfigures (a), (b), (c), and (d) present the respective fitted results of
the hybrid model, the neural network, their errors, and the training process,
respectively. From the results of these two experiments, we can conclude that
the hybrid model fit the data well, and the neural network model fit the
infection function well. Furthermore, the distribution of the errors
illustrate no obvious laws. Comparing Figure 5 with Figure 6, during the
training process, the process of estimating $\beta(S)$ had more vibrations.

##  6 Application to real data

In this section, we apply the proposed methods to predict the ground-truth
changes in historical COVID-19 infectivity. U.S. COVID-19 data were used,
collected from the official World Health Organization (WHO,
https://covid.cdc.gov/covid-data-tracker/#trends_dailycases) website. Time
periods, June 16 to October 31, 2021, and December 12, 2021, to February 11,
2022, were used for model fitting. These periods present the two most recent
case peaks in 2021–2022. During COVID-19 transmission, the dominant strains
were delta and omicron, respective to the two peaks [15]. To reduce the impact
of decreased case detection on weekends and the subsequent increases on
Mondays, a 7-day moving average was used. The data are available from the WHO
website.

(a) Fitted result of the hybrid model of Equation (8) during the first period.
(b) Fitted result of the hybrid model of Equation (8) during the second
period.

Figure 7: Fitted results of the model of Equation (8) with two datasets
collected from (a) June 16 to October 31, 2021, and (b) December 12, 2021, to
February 11, 2022, from U.S. COVID-19 case data.

To appropriately fit the data so that the infection function properties can be
studied, we introduce a modified SIR model:

$\left\\{\begin{aligned} \frac{d\hat{S}}{dt}&=-\beta\hat{S}I,\\\
\frac{dI}{dt}&=\beta\hat{S}I-(\gamma+d)I,\end{aligned}\right.$

where $\hat{S}=S\times 10^{-3}$; parameter $\gamma$ is the recovery rate, and
$d$ is the death rate. In the model, the $\hat{S}$ is a variable of the
infection function, $\beta\hat{S}I$, which shares its trend with $\beta SI$ of
the standard SIR model. We also introduce a sufficiently large constant, $M$,
such that for all times $t$, $\beta_{M}=M\beta$, $s=\hat{S}/M$, and $i=I/M$
follow $s(t)<1$ and $i(t)<1$ to produce the following model equation:

$\left\\{\begin{aligned} \frac{ds}{dt}&=-\beta_{M}si,\\\
\frac{di}{dt}&=\beta_{M}si-(\gamma+d)i.\end{aligned}\right.$ (8)

The initial values were $s_{0}=1-10^{-3}i_{0}$, $i_{0}=I_{0}/M$, and
$i(t)=I(t)/M$. For calculation, we used $M=7.6N\times 10^{-3}$, where $N$ is
the population of the U.S. During transmission, temperature changes and
behavioral changes affected the infection rates of the two key virus strains.
Infection rate $\beta$ is a time-related function, $\beta(t)$. Thus, in the
above model of Equation (8), we can assume that $\beta=\beta(t)$ is the
function that must be fitted by the neural network.

Figure 8: Fitted result of infection function $\beta(t)$ of the model of
Equation (8) using data collected from two time periods around the two peaks
of the COVID-19 case data.

Figures 7 (a) and (b) illustrate the model-fitted results of first and second
peak U.S. COVID-19 case data, and Figure 8 shows the fitted results of the
infection function of these two peaks. In the first period, the infectivity
first decreased and then became constant. In the second period, the
infectivity first increased and then became constant. However, after comparing
the fitted results of the two infection functions, we can conclude that the
change rate of the infection function in the second period was faster than
that of the first. Similarly, the second period had a higher infection rate
than the first. The decrease in the first period may have been caused by the
increased number of vaccinated people, and the increase in the second may have
been caused by the introduction omicron strain and the decrease in
temperature.

Through this practical application example, it was demonstrated that our
method can interpret otherwise hidden information in limited real-world data,
which not only increases application breadth but also greatly increases model
adaptability. The proposed method establishes a theoretical basis for further
research on infectious diseases and provides algorithmic support for the same.

##  7 Conclusion and discussion

In this work, by introducing the basic properties of an epidemic model and the
forward bifurcation of an epidemic system, we explained the vanishing gradient
conditions that result in the untrainability of the hybrid model. A novel loss
function that combines MSE loss and bifurcation items was then constructed to
overcome this issue. The existence conditions needed to ensure the
estimability of the neural network were also discovered. Based on the standard
fixed-step numerical method, we trained our hybrid model with generated data,
and the numerical results verify its accurate estimability. To our knowledge,
this is the first work to theoretically analyze and empirically demonstrate a
trainable nonlinear hybrid model that accurately predicts the behavior of
real-world behaviors.

In Section 5, the results of four numerical experiments were described to
demonstrate the applicability and efficacy of the hybrid model in estimating
periodical and Holling type infection functions. All experiments show that the
system of Equation (1) can fit the data very well, validating Theorem 2. We
also provided the necessary parameters to complete the experiments alongside a
specific implementation method and its published programming code, all of
which may be used in related hybrid-model applications. In Section 6, we
reported the application of our method using real COVID-19 data from the U.S.,
and the results showed that the proposed model indeed finds the hidden
information behind limited real data, which increases its applicability to
infectious disease models.

This work provides exciting new sights into the capabilities of hybrid neural
networks in handling nonlinear problems via nonlinear ordinary differential
equations. In a future work, we plan to focus on more complex but more widely
used stochastic and partial differential equations. Such follow-on studies
will extend the range of the application of hybrid models while enabling them
to solve complex econometric, computational fluid mechanic, and computational
chemistry problems. Furthermore, deep-learning models’ great dependency on
manually labeled training data may be mitigated. As the hybrid model is based
on many prior conditions, it can rely on smaller neural networks, which
reduces the need for training data.

## References

  * [1] ME Alexander and SM Moghadas.  Periodicity in an epidemic model with a generalized non-linear incidence.  Mathematical Biosciences, 189(1):75–96, 2004.
  * [2] Roy M Anderson and Robert M May.  Infectious diseases of humans: dynamics and control.  Oxford university press, 1992.
  * [3] Jingan Cui, Yonghong Sun, and Huaiping Zhu.  The impact of media on the control of infectious diseases.  Journal of dynamics and differential equations, 20(1):31–53, 2008\.
  * [4] Xiangjun Du, Aaron A King, Robert J Woods, and Mercedes Pascual.  Evolution-informed forecasting of seasonal influenza a (h3n2).  Science translational medicine, 9(413), 2017.
  * [5] Bing Guo, Dingkai Li, Congming Cheng, Zi-an Lü, and Youting Shen.  Simulation of biomass gasification with a hybrid neural network model.  Bioresource Technology, 76(2):77–83, 2001.
  * [6] Francis Begnaud Hildebrand.  Introduction to numerical analysis.  Courier Corporation, 1987.
  * [7] Tim Hill, Leorey Marquez, Marcus O’Connor, and William Remus.  Artificial neural network models for forecasting and decision making.  International journal of forecasting, 10(1):5–15, 1994.
  * [8] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366, 1989.
  * [9] Lars Hufnagel, Dirk Brockmann, and Theo Geisel.  Forecast and control of epidemics in a globalized world.  Proceedings of the National Academy of Sciences, 101(42):15124–15129, 2004.
  * [10] Clifford M Hurvich and Chih-Ling Tsai.  Regression and time series model selection in small samples.  Biometrika, 76(2):297–307, 1989.
  * [11] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [12] Kun Chang Lee, Ingoo Han, and Youngsig Kwon.  Hybrid neural network models for bankruptcy predictions.  Decision Support Systems, 18(1):63–72, 1996.
  * [13] Zhien Ma, Yicang Zhou, and Jianhong Wu.  Modeling and dynamics of infectious diseases, volume 11.  World Scientific, 2009.
  * [14] Maia Martcheva.  Methods for deriving necessary and sufficient conditions for backward bifurcation.  Journal of biological dynamics, 13(1):538–566, 2019.
  * [15] Ranjan K Mohapatra, Ruchi Tiwari, Ashish K Sarangi, Sanjay K Sharma, Rekha Khandia, G Saikumar, and Kuldeep Dhama.  Twin combination of omicron and delta variant triggering a tsunami wave of ever high surges in covid-19 cases: a challenging global threat with a special focus on indian sub-continent.  Journal of medical virology, 2022.
  * [16] Dimitris C Psichogios and Lyle H Ungar.  A hybrid neural network-first principles approach to process modeling.  AIChE Journal, 38(10):1499–1511, 1992.
  * [17] Navid Razmjooy, Fatima Rashid Sheykhahmad, and Noradin Ghadimi.  A hybrid neural network–world cup optimization algorithm for melanoma detection.  Open Medicine, 13(1):9–16, 2018.
  * [18] Yafeng Ren, Hao Fei, Xiaohui Liang, Donghong Ji, and Ming Cheng.  A hybrid neural network model for predicting kidney disease in hypertension patients based on electronic health records.  BMC medical informatics and decision making, 19(2):131–138, 2019\.
  * [19] Shigui Ruan and Wendi Wang.  Dynamical behavior of an epidemic model with a nonlinear incidence rate.  Journal of Differential Equations, 188(1):135–163, 2003.
  * [20] HL Smith.  Multiple stable subharmonics for a periodic epidemic model.  Journal of Mathematical Biology, 17(2):179–190, 1983.
  * [21] Dennis Wackerly, William Mendenhall, and Richard L Scheaffer.  Mathematical statistics with applications.  Cengage Learning, 2014.
  * [22] Wolfgang Walter.  Ordinary differential equations.  182\. Springer, 1998.
  * [23] Wendi Wang and Xiao-Qiang Zhao.  Threshold dynamics for compartmental epidemic models in periodic environments.  Journal of Dynamics and Differential Equations, 20(3):699–717, 2008\.
  * [24] Wendi Wang and Xiao-Qiang Zhao.  Basic reproduction numbers for reaction-diffusion epidemic models.  SIAM Journal on Applied Dynamical Systems, 11(4):1652–1673, 2012\.
  * [25] Hong Wu, Yan Xing, Weihong Ge, Xiaoquan Liu, Jianjun Zou, Changjiang Zhou, and Jun Liao.  Drug-drug interaction extraction via hybrid neural networks on biomedical literature.  Journal of biomedical informatics, 106:103432, 2020.
  * [26] Dongmei Xiao and Shigui Ruan.  Global analysis of an epidemic model with nonmonotone incidence rate.  Mathematical biosciences, 208(2):419–429, 2007.
  * [27] Kiyoshi Yamaoka, Terumichi Nakagawa, and Toyozo Uno.  Application of akaike’s information criterion (aic) in the evaluation of linear pharmacokinetic equations.  Journal of pharmacokinetics and biopharmaceutics, 6(2):165–175, 1978.
  * [28] Toshie Yamashita, Keizo Yamashita, and Ryotaro Kamimura.  A stepwise aic method for variable selection in linear regression.  Communications in Statistics—Theory and Methods, 36(13):2395–2403, 2007.
  * [29] Fang Zhang and Xiao-Qiang Zhao.  A periodic epidemic model in a patchy environment.  Journal of Mathematical Analysis and Applications, 325(1):496–516, 2007.
  * [30] Liu Ziyin, Tilman Hartwig, and Masahito Ueda.  Neural networks fail to learn periodic functions and how to fix it.  Advances in Neural Information Processing Systems, 33, 2020.


* (82) R. R. Landim, “On the Laplace equation in d-dimension,” [arXiv:0901.3488 [math-ph]].