=Paper=
{{Paper
|id=Vol-3831/paper13
|storemode=property
|title=Mechanistic Causal Models for Explainable AI in Medicine: Coupling Respiratory and Immunological Systems for In Silico Medicine Simulations
|pdfUrl=https://ceur-ws.org/Vol-3831/paper13.pdf
|volume=Vol-3831
|authors=José Ignacio Lorenzo,María Dolores Corbacho,Fernando Corbacho
|dblpUrl=https://dblp.org/rec/conf/explimed/LorenzoCC24
}}
==Mechanistic Causal Models for Explainable AI in Medicine: Coupling Respiratory and Immunological Systems for In Silico Medicine Simulations==
https://ceur-ws.org/Vol-3831/paper13.pdf
Mechanistic Causal Models for Explainable AI in
Medicine: Coupling Respiratory and Immunological
Systems for In Silico Medicine Simulations⋆
José Ignacio Lorenzo1 , María Dolores Corbacho2 and Fernando Corbacho3,∗
1
Universidad Autónoma de Madrid, Ctra. Colmenar Viejo Km 15, Madrid, 28049, Spain
2
Hospital Ribera Povisa, Rúa de Salamanca 5, Vigo, 36211, Spain
3
Cognodata R&D, 𝑃.𝑜 de la Castellana 135, Madrid, 28046, Spain
Abstract
Explainable AI for the medical domain is a critical necessity. We propose explainable mechanistic causal
artificial intelligence models that incorporate known physiological, anatomical, physical, chemical and
control principles while, at the same time, allowing data driven parameter estimation to personalize to
specific patient´s personal characteristics. This type of hybrid systems medicine approach, also incorpo-
rating control theory principles, not only allows personalized simulations but also causal explanations
identifying and explaining irregular clinical states. While pure data-driven systems can create highly
accurate models, there are concerns about their fairness, opacity and lack of explainability. Therefore, we
promote hybrid approaches that enhance transparency, accountability, trustworthiness and explainability.
In this paper, we focus on mechanistic causal models of the respiratory and the immunological systems
incorporating also control theory principles to explain the dynamics a particular pathology, namely, the
cytokine storm. Yet, the methodology here proposed, is applicable to any medical domain, as well as, to
the coupling of different medical systems. Thus, helping to bridge the gap towards more explainable
systems medicine.
Keywords
Simulation, Homeostasis, Feedback and Feedforward control, Stability, Clinical state, Personalized
Systems Medicine.
1. Introduction
We propose explainable mechanistic causal artificial intelligence models that incorporate known
physiological, anatomical, physical, chemical and control principles. Mechanistic models are
mainly based on the causal understanding of biological entities (e.g. proteins, cells, organs)
and their dynamic interactions [1]. Basically, these models represent existing knowledge about
biological systems in a form that can generate predictions on the behavior of the system. They
also allow to explore biomedical simulations for personalized systems medicine. Thus, mech-
anistics models are white-box, interpretable, models where any dynamic behavior, including
the pathological behaviors, can be explained, as shown in this article. As Assaad et al. [2]
EXPLIMED - First Workshop on Explainable Artificial Intelligence for the medical domain - 19-20 October 2024, Santiago
de Compostela, Spain.
∗
Corresponding author.
Envelope-Open josei.lorenzo@estudiante.uam.es (J. I. Lorenzo); dcorbacho@povisa.es (M. D. Corbacho);
fernando.corbacho@cognodata.com (F. Corbacho)
© 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR
ceur-ws.org
Workshop ISSN 1613-0073
Proceedings
emphasize, causality is indeed crucial for explanatory purposes since, thanks to it, an effect can
be explained by its causes, regardless of the correlations it may have with other variables.
While pure data-driven systems can create highly accurate models, there are concerns about
their fairness, opacity and lack of explainability. Therefore, we promote hybrid approaches
that enhance transparency, accountability, trustworthiness and explainability. This type of
hybrid systems medicine approach, also incorporating control theory principles, not only allows
personalized simulations but also causal explanations identifying and explaining irregular
clinical states. In this paper, we focus on mechanistic causal models of the respiratory and the
immunological systems incorporating also control theory principles to explain the dynamics a
particular pathology, namely, the cytokine storm.
In this paper, we focus on two vital systems: the human respiratory system and the immune
system. Using simple equations to better identify the variables that, possibly, cause irregularities
and more advanced control theory techniques, we investigate how these systems interact with
each other and respond to extreme conditions, such as cytokine storms. The human respiratory
system has been studied and attempted to reproduce e.g. [3]. It’s crucial for maintaining normal
oxygen and carbon dioxide pressures in humans. Its main function is to adapt pulmonary
ventilation to the metabolic needs of consumption and production of both gases. Despite
variations in oxygen requirements and carbon dioxide elimination, the pressure of these elements
in the lungs is maintained within narrow ranges thanks to a complex regulation of pulmonary
ventilation [4]. So it will start from the pressure equation to develop the different parts that
make up the respiratory system, and it will see what changes occur when the average oxygen
concentration needs to be changed.
The brainstem includes a Central Pattern Generator (CPG) that regulates breathing [5]. The
respiratory biological neural network changes the respiratory pattern, which can be controlled
by the frequency and amplitude to maintain physiological levels of oxygen [6]. In fact, there
are many mixtures of them that can achieve the same physiological levels [7]. So the question
is: can a simulation of the respiratory system controlled by the amplitude and frequency be
carried out to predict the evolution of a patient? Other papers have described more biologically
plausible neural networks for Central Pattern Generators (CPGs) and MPGs in other systems [8],
[9] and in the respiratory system [10]. Yet, in this paper, for simplicity reasons, we have opted
for a more phenomenological dynamical description by a sort of ”sinusoidal” function. For this
reason, the equations presented in [11] have been used, yet these do not take into account the
muscular activity in the lungs, so an equation has been added to better study the predictive
simulation. In addition, several concepts of control theory will also be used. Specifically, a PID
feedback controller and inverse model for feedforward control will be used. All the differential
equations have been solved numerically for each instant of time using the Runge-Kutta method
of order 4.
Additionally, a simulation of the immune system is also carried out, highlighting its simplicity,
as well as, its adaptability against pathogenic agents. Once this is done, it is coupled with the
simulated respiratory system and their interactions are analyzed to maintain the homeostatic
regime. Finally, the effect that a cytokine storm produces on both systems is simulated, re-
sulting in a disproportionate immune response that causes the respiratory system to fail. This
methodology, not only seeks to better understand and explain the diseases, but, eventually, also
to develop more effective therapeutic strategies.
2. Respiratory system equations and parameters
The equations and values of the variables are inspired by [11] with some adaptations. Also, an
equation that describes muscle activity has been included. The human respiratory system’s
main objective is to maintain a normal oxygen concentration, which is achieved among other
things, through the CPG’s amplitude and period parameters. The CPG is a biological neural
network located in the brainstem responsible for generating the basic respiratory rhythm, and
coordinating the contraction and relaxation of respiratory muscles, such as the diaphragm and
intercostal muscles. The 𝑝(𝑡) function in [11] refers to lung pressure, but when introducing
muscle action in this article, it represents the CPG output control signal that acts on the
respiratory muscles, as it is depicted in the following equations. The pressure 𝑝(𝑡) is defined as
follows:
𝐴 2𝜋𝑅
( 𝑇 𝑠𝑖𝑛 8𝜋𝑡 − 𝐸𝑐𝑜𝑠 8𝜋𝑡 ) + 𝐸𝐴 + 𝐸𝐵 if 𝑡 ∈ (0, 3𝑇 ]
𝑝(𝑡) = { 𝐴2 2𝜋𝑅 3𝑇 3𝑇 2 8 (1)
(
2 𝑇
𝑠𝑖𝑛( 5𝑇 + 5 )𝜋 − 𝐸𝑐𝑜𝑠( 5𝑇 + 25 )𝜋) + 𝐸𝐴
8𝑡 2 8𝑡
2
+ 𝐸𝐵 3𝑇
if 𝑡 ∈ ( 8 , 𝑇 ]
where 𝐵 = 0.2 𝑙 is the minimum lung volume, 𝐸 = 2.5 𝑠 −1 is the elasticity of the lungs, 𝑅 = 1
𝑚𝑚𝐻 𝑔 ⋅ 𝑠 ⋅ 𝑙 −1 is the resistance to airflow, 𝑡 is time, 𝑇 is the period in the instant of time 𝑡 and
𝐴 the amplitude. Furthermore, it is considered that the respiration is divided in two phases:
inspiration occurs when 𝑡 ∈ (0, 3𝑇 8
) and expiration occurs when 𝑡 ≥ 3𝑇 8
. The muscular activity
𝑚(𝑡) in the instant of time 𝑡 is:
𝑚′ (𝑡) = −𝑘1 𝑚(𝑡) + 𝑘2 𝑝(𝑡) (2)
where 𝑘1 = 2 is the recoil rate constant of muscle, 𝑘2 = 1 is the activation rate constant of
muscle. The lung volume 𝑣(𝑡) in the instant of time 𝑡 is:
1
𝑣 ′ (𝑡) =
[−𝐸𝑣(𝑡) + 𝑚(𝑡)] (3)
𝑅
Finally, the oxygen concentration 𝑜(𝑡) at time 𝑡 is:
𝐴−𝐷𝑆 ′
−𝐷𝑜(𝑡) + 𝐼 ( )𝑣 (𝑡) if 𝑡 ∈ (0, 3𝑇 ]
𝑜 ′ (𝑡) = { 𝐴
3𝑇
8 (4)
−𝐷𝑜(𝑡) if 𝑡 ∈ ( 8 , 𝑇 ]
where 𝐷 = 0.75 𝑠 −1 is the diffusion rate of 𝑂2 , 𝑁 = 0.15 𝑙 is the death space or volume of
the airways that does not participate in gas exchange, 𝐼 = 2.25 𝑙 −1 is the rate of increase in
oxygen concentration and 𝑜𝑎𝑣 is the average value of oxygen concentration. Figure 1 depicts the
main components of the respiratory system, namely, the CPG, the lungs as compartments, the
muscles as springs, the oxygen concentration chemoreceptors, as well as, the control structures
to be introduced in later sections.
3. PID Controller & Simulations
To begin with, a feedback controller will be added to study different respiratory control mecha-
nisms. Similarly to [10], who inlude a closed-loop respiratory system model of the brainstem
Figure 1: Respiratory system with adaptative control. Where 𝑝(𝑡) represents the pressure, 𝑣(𝑡) the
pulmonary volume, 𝑚(𝑡) the muscular activity and 𝑜(𝑡) the oxygen concentration.
respiratory network that controls the pulmonary system, and a subsystem that represents
lung biomechanics, and gas exchange and transport (𝑂2 and 𝐶𝑂2 ). The biological pulmonary
subsystem provides two types of feedback to the neural subsystem: a mechanical one coming
from the lung stretch receptors, and another chemical one coming from the chemoreceptors.
Yet, [10] focus more on simulating respiratory neurons that interact within the Botzinger and
pre-Botzinger complexes, as well as the retrotrapezoidal parafacial respiratory nucleus/group
(RTN/pFRG) representing the central chemoreception module targeted by the chemical. Also,
[12] presents a system combing clinical evidence and expert knowledge within a physiological
closed-loop control structure for mechanical ventilation.
Thus, simulation of the human respiratory system is a powerful tool for research and de-
velopment of respiratory therapies. A starting point for these simulations is the inclusion
of a Proportional-Integral-Derivative (PID) controller, which allows the system’s response to
changing conditions in the environment to be dynamically adjusted. A PID controller is a simple
control mechanism that, through a feedback loop, allows regulating pressure, muscle activity,
lung volume and oxygen concentration. The PID controller calculates the difference between
the current value for the variable versus the desired value for the same variable. The PID control
algorithm consists of three different components: proportional, integral, and derivative.
1. Proportional: = 𝐾𝑝 𝑒(𝑡), depends on the current error.
𝑡
2. Integral: = 𝐾𝑖 ∫0 𝑒(𝜏 )𝑑𝜏, depends on past errors.
𝑑𝑒(𝑡)
3. Derivative: = 𝐾𝑑 𝑑𝑡 , is a prediction of future errors.
where 𝑒(𝑡) it is a function that represents the error between the value calculated at a moment
in time and the target value that you want to achieve.
Previously, there are articles that apply a PID controller to the respiratory system such as
[3]. In this paper, the system modeled focuses on the average oxygen concentration, 𝑜𝑎𝑣 . At
∗ , and the controller attempts to achieve it by modifying either
each instant, there is a target 𝑜𝑎𝑣
the amplitude 𝐴, or the period 𝑇, or both. This change in the variables is made through a
conversion of the PID controller. If the controller only modifies 𝑇, then the change corresponds
to 𝑇 = 𝑇 − 𝑘 where k is the signal from the PID controller. In case that it only modifies 𝐴, the
change corresponds to 𝐴 = 𝐴 + 𝑘. However, if both are modified, then the change corresponds
to 𝐴 = 𝐴 + 0.1𝑘, 𝑇 = 𝑇 − 𝑘. Once the change is made, pressure, muscle activity, lung volume
and oxygen concentration correspondingly change.
Several simulation tests have been performed to validate the results. One of them is shown
below in Figure 2. During this simulation, 𝑜𝑎𝑣∗ is changed 9 times, with a duration of 10 seconds
for each change, starting from an initial 𝑜𝑎𝑣 ∗ = 0.08 mg/l, which in total makes the whole
simulation process last 100 seconds, as can be seen in Figure 2 with the blue line. On the other
hand, the orange line refers to the values of 𝑜𝑎𝑣 that are calculated by the controller. In this
specific simulation, the controller (only) modifies the amplitude 𝐴 so that 𝑜𝑎𝑣 can reach the
∗ . Table 1 includes 9 different changes in 𝑜 ∗ , one in each row. A change occurs
target value 𝑜𝑎𝑣 𝑎𝑣
every 10 seconds as displayed in Figure 2, and each change has a different magnitude. For each
change, each row indicates the size (magnitude) of the change, the number of iteractions, and
the time in seconds, that it took the PID controller to achieve that particular reference value of
∗.
𝑜𝑎𝑣
∗
Figure 2: The PID controller regulates 𝑜𝑎𝑣 (orange line) to approach 𝑜𝑎𝑣 (blue line) by varying A.
4. Feedforward Control & Inverse model
A PID controller is a control mechanism that, through a feedback loop, allows to calculate
the difference between the current value of the variable versus the desired value for the same
variable. Yet, it always acts a posteriori, that is, there must be a change in the real value that
you want to approximate for it to start acting. Most articles discuss this type of closed-loop
feedback control of the human respiratory system as in [3], [4], [10], [12]. On the other hand, a
feedforward controller can anticipate this change and act before the change takes place and,
thus, reduce the number of interactions needed by the controller. In this case, a feedforward
inverse model is designed to learn through a training set that approximates the objective faster
than the PID controller. Hence, the feedfoward inverse model allows for explainability of
anticipatory responses experimentally observed in the nervous system.
The inverse model is essentially a function that maps the inputs to the outputs. During the
training phase, the inverse model is exposed to known input and output data, which come
Table 1
On the left, the results of the PID Controller when changing only 𝐴. On the right, the results of the
Feedforward Inverse model when changing only 𝐴. Ten tests have been run for each specific change
∗
of target mean oxygen concentration 𝑜𝑎𝑣 . The column named change indicate when the change takes
place. The column named size refers to the magnitude of that specific change. We present the mean
and standard deviation resulting from the ten tests run for each specific case.
change size iterations time (s) change size iterations time (s)
mean std mean std mean std mean std
10 s 0.19 18.3 3.1 1.851 0.26 10 s 0.19 3.8 0.75 0.411 0.07
20 s -0.19 20 2.57 2.011 0.37 20 s -0.19 4.3 1 0.521 0.2
30 s 0.13 15.9 1.92 1.71 0.23 30 s 0.13 3.2 0.75 0.369 0.12
40 s 0.16 16.4 1.74 1.744 0.21 40 s 0.16 4.7 2.15 0.513 0.25
50 s -0.2 20.1 2.95 2.259 0.46 50 s -0.2 3.9 0.83 0.433 0.11
60 s -0.02 2.5 2.16 0.253 0.22 60 s -0.02 1.8 0.75 0.191 0.09
70 s 0.1 13.5 2.16 1.496 0.28 70 s 0.1 3.8 1.94 0.457 0.19
80 s -0.04 6.9 2.12 0.723 0.19 80 s -0.04 1.7 0.9 0.19 0.09
90 s -0.14 16.1 2.07 1.754 0.17 90 s -0.14 3 0.63 0.398 0.14
mean 14.4 2.31 1.533 0.27 mean 3.4 1.08 0.387 0.14
from control tests performed on the actuator. The inverse model adjusts its internal weights
to minimize the error between the predicted outputs and the actual outputs of the actuator.
Once trained, the inverse model can predict the actuator outputs for new inputs, even if they
have not been seen before. After training, the learned inverse model can be used to control
the actuator. When the inverse model only modifies the variable 𝑇 of section 2, the change
corresponds to 𝑇 = 𝑇 − 10𝑘, where k is the signal of inverse model. In case it only modifies 𝐴,
the change corresponds to 𝐴 = 𝐴 + 𝑘. However, if both are modified by the inverse model, then
𝐴 = 𝐴 + 0.75𝑘, 𝑇 = 𝑇 − 7.5𝑘.
4.1. Feedback Error Learning
The feedback and feedforward controllers applied to human motor systems are needed. In this
paper, as in [13], first a feedback controller has been used that approximates the target value of
average oxygen concentration, and then the error is calculated to better anticipate the value
that the signal will give with a feedforward controller trained with feedback error learning. The
equations that describe the process are:
𝑓𝑓 𝑓𝑏
𝑥𝑡+1 = 𝑓 (𝑥𝑡 , 𝑢𝑡 ) = 𝑓 (𝑥𝑡 , 𝑢𝑡 + 𝑢𝑡 ) (5)
𝑓𝑏 𝑓𝑓
where 𝑢𝑡 is a signal from the PID controller (feedback), and 𝑢𝑡 is the signal from the
feedforward controller. 𝑓 corresponds to the plant dynamics and 𝑔, below, corresponds to the
inverse model parametrized by w.
𝑓𝑓 ∗ , 𝑥 ; w)
𝑢𝑡 = 𝑔(𝑥𝑡+1 𝑡 (6)
In Figure 1, it can be seen that the inverse model receives a current value and a target value
that it has to achieve. This creates a signal that affects the amplitude 𝐴, or period 𝑇, or both.
This, in turn, changes the CPG control signal 𝑝(𝑡) on the muscle 𝑚(𝑡), which also affects the
lung volume 𝑣(𝑡), and the oxygen concentration 𝑜(𝑡). Once the process is complete, the average
oxygen concentration changes and the inverse model learns about the error made to better
approximate the target values.
4.2. Linear Inverse Model & Parameters Update
At the beginning of training, the internal weights of the inverse model are randomly initialized.
These weights, w, are the parameters that the system will adjust to learn the inverse model.
During each training iteration, input data is provided to the inverse model, that performs
forward propagation, computing the control outputs using the current weights 𝑤𝑡𝑖 . In this
section we will describe a linear inverse model:
𝑓𝑓
𝑢𝑡 = 𝑤𝑡1 ⋅ 𝑜𝑎𝑣 + 𝑤𝑡2 ⋅ 𝑜𝑎𝑣
∗ (7)
where 𝑜𝑎𝑣 is the previous value of average oxygen concentration, and 𝑜𝑎𝑣∗ is the value of the
target average oxygen concentration. The learning algorithm calculates the partial derivatives
of the error with respect to each weight; these derivatives indicate how to change the weights
to reduce the error. The error signal in feedback error learning corresponds to the output of
the feedback controller [13]. Then, the weights are updated using the gradient descent rule;
basically the weights adjust in the direction that reduces the error, namely,
𝑑𝑔𝑡 𝑓𝑏
∇𝑤𝑡𝑖 = 𝜖 𝑖
𝑢𝑡 (8)
𝑑𝑤𝑡
As can be seen in TABLE 1 (right), with this inverse model and feedforward control, it is
possible to drastically reduce the number of iterations and, consequently, the time it takes to
reach the objective, compared to the results obtained in TABLE 1 (left). A nonlinear Inverse
model implemented by a multilayer neural network trained by backpropagation has also been
developed. It is not included due to space limitations and, also, due to the fact that it does not
improve the results obtained by the linear, simpler, better for explainability linear model.
5. Immune equations and parameters
Cytokine storm (section 5.1) is a life-threatening inflammatory response characterized by hy-
peractivation of the immune system and can be caused by various therapies, autoimmune
conditions, or pathogens. Several mathematical models have been developed to try to describe
the dynamics of cytokine storms. In [14], cytokines are divided into 7 categories. They use
a model with data from mice to describe the interactions of cytokines with each other. In
[15], they experimented with six volunteers who had a severe inflammatory response during
the clinical trial, which caused a cytokine storm that could be studied using a set of ordinary
differential equations. [16] applies a simple two-component differential equation model for
pro- and anti-inflammatory responses and detailed mathematical analysis to identify specific
responses to cytokine storms. [17] relates the cytokine storm as a cause of SARS-CoV. For this
purpose, a system of fifteen ordinary differential equations is presented that models the immune
response to SARS-CoV-2 that infects immune cells, which can trigger a cytokine storm.
However, [18] has been considered, for explainability reasons, most appropriate to combine
with the simulated respiratory system. In this model, the five most important variables are:
susceptible cells 𝑆(𝑡), infected cells 𝐼 (𝑡), viral particles 𝑉 (𝑡), immune cells 𝑥(𝑡) and the cytokines
𝑦(𝑡). Susceptible cells have a turnover that is reflected in 𝑆𝑖𝑛 − 𝑘𝑆 𝑆(𝑡) and can be infected with a
constant 𝛽𝑖𝑠 , represented by 𝛽𝑉 (𝑡)𝑆(𝑡). Infected cells are eliminated at a rate 𝑘𝐼 naturally 𝑘𝐼 𝐼 (𝑡),
and at a rate 𝛾 by immune cells 𝛾 𝑥(𝑡)𝐼 (𝑡). Viral particles are produced by infected cells 𝑣𝑖𝑛 𝐼 (𝑡)
and eliminated 𝑘𝑣 𝑉 (𝑡). These equations have been reproduced in Figures 3 and it has been
verified that the results are consistent:
𝑑𝑆
= 𝑆𝑖𝑛 − 𝑘𝑠 𝑆(𝑡) − 𝛽𝑉 (𝑡)𝑆(𝑡) (9)
𝑑𝑡
𝑑𝐼
= 𝛽𝑉 (𝑡)𝑆(𝑡) − 𝑘𝐼 𝐼 (𝑡) − 𝛾 𝑥(𝑡)𝐼 (𝑡) (10)
𝑑𝑡
𝑑𝑉
= 𝑣𝑖𝑛 𝐼 (𝑡) − 𝑘𝑣 𝑉 (𝑡) (11)
𝑑𝑡
In the previous system, the interaction between immune cells 𝑥(𝑡) and cytokines 𝑦(𝑡) is
not described. The normal transformation of immune cells 𝑥(𝑡) is given by 𝑥𝑖𝑛 − 𝑘1𝑖𝑠 𝑥(𝑡), the
production induced by infected cells 𝛾1 𝑥(𝑡)𝐼 (𝑡) and the production of additional immune cells
𝑥(𝑡)
by 𝑏1 𝑐 +𝑥(𝑡) (𝑥(𝑡) − 𝑚𝑖𝑠 )(𝑦1 − 𝑦(𝑡))(𝑦(𝑡) − 𝑦2 ). The normal transformation of cytokines is given by
1
𝑎 𝑦(𝑡)𝑥(𝑡)
𝑦𝑖𝑛 − 𝑘2𝑖𝑠 𝑦(𝑡) and the production stimulated by immune cells 𝑏2 𝑎 1(𝑐 +𝑥(𝑡)) . Thus, the followings
2 2
equations have been reproduced in Figures 3 :
𝑑𝑥 𝑥(𝑡)
= 𝑥𝑖𝑛 − 𝑘1𝑖𝑠 𝑥(𝑡) + 𝛾1𝑖𝑠 𝑥(𝑡)𝐼 (𝑡) + 𝑏1 (𝑥(𝑡) − 𝑚𝑖𝑠 )(𝑦1 − 𝑦(𝑡))(𝑦(𝑡) − 𝑦2 )) (12)
𝑑𝑡 𝑐1 + 𝑥(𝑡)
𝑑𝑦 𝑎1 𝑦(𝑡)𝑥(𝑡)
= 𝑦𝑖𝑛 − 𝑘2𝑖𝑠 𝑦(𝑡) + 𝑏2 (13)
𝑑𝑡 𝑎2 (𝑐2 + 𝑥(𝑡))
5.1. Cytokine storm
A cytokine storm occurs when there is a positive feedback between cytokines and immune cells.
Cytokines direct immune cells to the site of infection and stimulate the production of more
cytokines, however, IL-6 and IL-10 interleukins break the natural transition from inflammation
to recovery. Cytokine storm can cause significant damage to the epithelial cells of the lung
alveoli. It has been established that the value of viral particles is 𝑉 (500) = 0.8, which produces
a cytokine storm. Once this is done, we can observe and explain (very important for explainable
AI in Medicine) the effect it has on the simulated respiratory system. In a cytokine storm,
certain proinflammatory interleukins act to damage airway epithelial cells and inhibit oxygen
uptake in the lungs, so the following variables of the respiratory system seen in Section 2 have
been modified: the diffusion rate of oxygen concentration 𝐷 and the rate of increasing oxygen
concentration 𝐼, that affect the amplitude, period, or both, to achieve the target 𝑜𝑎𝑣 .
𝑦(𝑡𝑖+1 ) − 𝑦(𝑡𝑖 )
𝐷(𝑡𝑖+1 ) = 𝐷(𝑡𝑖 ) + 0.1 (14)
3
Figure 3: 𝑆 is the volume of susceptible cells measured as a percentage/ml in seconds, the number of
cells decreases significantly because they are infected. 𝐼 is the volume of infected cells measured as a
percentage/ml in seconds, the number increases significantly because the viral particles infect susceptible
cells. 𝑥 is the volume of immunes cells measured as a percentage/ml in seconds, the number increases
significantly when the cytokine storm ocurrs and decreases when the cytokine storm is finishing. 𝑦 is
the Cytokines storm has been simulated at 𝑡 = 500𝑠, where cytokines are measured in UI because in
[18] there is not a specific unit.
𝐼 (𝑡𝑖+1 ) = 𝐼 (𝑡𝑖 ) − 0.1(𝑦(𝑡𝑖+1 ) − 𝑦(𝑡𝑖 )) (15)
Figure 4: Up to the left, the amplitude 𝐴 reaches very high values in order to maintain 𝑜𝑎𝑣 during a
cytokine storm. The image in the upper right side is where pressure 𝑝(𝑡) is in common values before
cytokine storm occurs. After, at bottom left, the increase in amplitude 𝐴 causes a decontrol in pressure
𝑝(𝑡) during the cytokine storm. Finally, at bottom right, pressure 𝑝(𝑡) returns to common values after
cytokine storm.
As can be seen from the results in Figure 4, during the cytokine storm, the amplitude 𝐴 needs
to reach very high values to maintain the target oxygen concentration. This causes variables
such as pressure (Figure 4) to also increase significantly to values that are difficult, or even
impossible, to reach by the human physiological system, instead of common values (Figures
4). Therefore, complications and respiratory distress, derived from the action of cytokines, are
obtained by these mechanistic models, and most importantly, can be explained and provide
feasible values that can be supported by each specific patient. In this case, cytokines that are
stimulated in excess by the immune cells, create an irregular pathological clinical state. This
pathological state can be dealt with means that inhibit the cytokines or/and inhibit the immune
cells pathways. Current innovative mHealth technology that allows for monitoring of, for
instance, blood oxygen saturation [19] provides an exciting bridge for the applicability of the
models developed in this paper in real life medical applications.
6. Conclusions
In summary, the variables that have been simulated for the respiratory system are: CPG ampli-
tude and/or period, pressure, muscle activity, lung volume and average oxygen concentration.
Additionally, for the immune system, the variables that have been simulated are: susceptible
cells, infected cells, viral particles, immune cells and cytokines. This paper only shows the
results by modifying the amplitude, yet, the best results are obtained when the amplitude and
the period of the CPG are modified simultaneously. The effect of a cytokine storm pathology
on these systems has been analyzed, but this can be expanded to simulate more diseases, and
in turn, analyze which parameters would recover a patient’s homeostatic regime. Therefore,
complications and respiratory distress, derived from the action of cytokines, are obtained by
these mechanistic models, and most importantly, can be explained and provide feasible values
that can be supported by each specific patient. Current innovative mHealth technology that
allows for monitoring of, for instance, blood oxigen saturation [19] provides an exciting bridge
for the applicability of the models developed in this paper in real life medical applications. This
methodology, not only seeks to better understand and explain the diseases, but, eventually, also
to develop more effective therapeutic strategies. Thus, this is a first step in simulating the effect
of different pathologies, and it will be the basis to also simulate the effect of different (dynamical)
treatment regimes. It would also be quite interesting to be able to include the cardiovascular
system within the overall system described in this paper. The works that seem most appropriate
are [20], [21], [22], but they have a large number of equations and variables, so at the moment,
it is not easy to implement them, keeping the explainability, within the system that has already
been created.
Acknowledgment
Cognodata has received support from red.es (NextGenerationEU funds) for the execution
of the project titled: “In-silico Medicine con sistemas generativos de inteligencia artificial y
schema-based machine learning: pilotos en enfermedades infecciosas” and expedient number:
2021/C005/00145600. The clinical part of this research, is conducted as part of the doctoral work
carried at the epidemiology department of the University of Santiago de Compostela by MDC.
We also thank two anonymous referees for very helpful suggestions.
References
[1] F. Franssen, P. Alter, N. Bar, B. Benedikter, S. Iurato, D. Maier, M. Maxheim, F. Rössler,
M. Spruit, C. Vogelmeier, E. Wouters, B. Schmeck, Personalized medicine for patients with
copd: where are we?, International Journal of Chronic Obstructive Pulmonary Disease 14
(2019) 1465–1484. doi:10.2147/COPD.S175706 .
[2] C. K. Assaad, E. Devijver, E. Gaussier, Survey and evaluation of causal discovery methods
for time series, J. Artif. Int. Res. 73 (2022). URL: https://doi.org/10.1613/jair.1.13428. doi:10.
1613/jair.1.13428 .
[3] A. Ben‐Tal, J. Smith, Control of breathing: Two types of delays studied in an integrated
model of the respiratory system, Respiratory physiology & neurobiology 170 (2010) 103–12.
doi:10.1016/j.resp.2009.10.008 .
[4] T. Miyamoto, System physiology of respiratory control in man, The Journal of Physical
Fitness and Sports Medicine 5 (2016) 329–337. doi:10.7600/jpfsm.5.329 .
[5] Y. Molkov, N. Shevtsova, C. Park, A. Ben-Tal, J. Smith, J. Rubin, I. A. Rybak, Computational
models of the neural control of breathing, Wiley interdisciplinary reviews. Systems biology
and medicine 9 (2017). doi:10.1002/wsbm.1371 .
[6] A. Ben‐Tal, M. Tawhai, Integrative approaches for modeling regulation and function of
the respiratory system, Wiley Interdisciplinary Reviews: Systems Biology and Medicine 5
(2013). doi:10.1002/wsbm.1244 .
[7] M. Tipton, A. Harper, J. Paton, J. Costello, The human ventilatory response to stress: Rate
or depth?, The Journal of Physiology 595 (2017) 5729–5752. doi:10.1113/JP274596 .
[8] F. Corbacho, K. Nishikawa, A. Weerasuriya, J.-S. Liaw, M. Arbib, Schema-based learning
of adaptable and flexible prey-catching in anurans i. the basic architecture, Biological
cybernetics 93 (2005) 391–409. doi:10.1007/s00422- 005- 0013- 0 .
[9] R. Huerta, M. Sánchez-Montañés, F. Corbacho, J. Siguenza, A central pattern generator
to control a pyloric-based system, Biological cybernetics 82 (2000) 85–94. doi:10.1007/
PL00007963 .
[10] Y. Molkov, N. Shevtsova, C. Park, A. Ben-Tal, J. Smith, J. Rubin, I. A. Rybak, A closed-loop
model of the respiratory system: Focus on hypercapnia and active expiration, PLoS ONE
9 (2014) e109894. doi:10.1371/journal.pone.0109894 .
[11] F. Zaidi, A. Ben-Tal, M. Roberts, Is our breathing optimal? solving a piecewise lin-
ear model with constraints, Journal of Mathematical Biology 83 (2021). doi:10.1007/
s00285- 021- 01661- 8 .
[12] P. von Platen, P. Pickerodt, M. Russ, M. Taher, L. Hinken, W. Braun, R. Köbrich,
A. Pomprapa, R. Francis, S. Leonhardt, M. Walter, Solve: a closed-loop system fo-
cused on protective mechanical ventilation, BioMedical Engineering OnLine 22 (2023).
doi:10.1186/s12938- 023- 01111- 0 .
[13] D. Wolpert, M. Kawato, Multiple paired forward and inverse models for motor control,
Neural Networks 11 (1998) 1317–1329.
[14] M. Waito, S. Walsh, A. Rasiuk, B. Bridle, A. Willms, A mathematical model of cytokine dy-
namics during a cytokine storm, Mathematical and Computational Approaches in Advanc-
ing Modern Science and Engineering (2016) 331–339. doi:10.1007/978- 3- 319- 30379- 6_
31 .
[15] H. Yiu, A. Graham, R. Stengel, Dynamics of a cytokine storm, PloS one 7 (2012) e45027.
doi:10.1371/journal.pone.0045027 .
[16] W. Zhang, S. Jang, C. Jonsson, L. Allen, Models of cytokine dynamics in the inflammatory
response of viral zoonotic infectious diseases, Mathematical medicine and biology : a
journal of the IMA 36 (2018). doi:10.1093/imammb/dqy009 .
[17] R. Reis, A. Pigozzo, C. Bonin, B. Quintela, L. Pompei, A. Vieira, L. Silva, M. Xavier, R. Santos,
M. Lobosco, A validated mathematical model of the cytokine release syndrome in severe
covid-19, Frontiers in Molecular Biosciences 8 (2021) 39423. doi:10.3389/fmolb.2021.
639423 .
[18] I. Kareva, F. Berezovskaya, G. Karev, Mathematical model of a cytokine storm, bioRxiv :
the preprint server for biology (2022). doi:10.1101/2022.02.15.480585 .
[19] G. Casalino, G. Castellano, G. Zaza, A mhealth solution for contact-less self-monitoring of
blood oxygen saturation, IEEE Symposium on Computers and Communications (ISCC)
(2020) 1–7. doi:10.1109/ISCC50000.2020.9219718 .
[20] A. Albanese, L. Cheng, M. Ursino, N. Chbat, An integrated mathematical model of the
human cardiopulmonary system: Model development, American Journal of Physiology -
Heart and Circulatory Physiology 310 (2016) ajpheart.00230.2014. doi:10.1152/ajpheart.
00230.2014 .
[21] E. Magosso, M. Ursino, Cardiovascular response to dynamic aerobic exercise: A
methematical model, Medical & biological engineering & computing 40 (2002) 660–74.
doi:10.1007/BF02345305 .
[22] C. Riley, A mathematical model of cardiovascular and respiratory dynamics in humans
with transposition of the great arteries, Rose-Hulman Undergraduate Mathematics Journal
18 (2017). URL: https://scholar.rose-hulman.edu/rhumj/vol18/iss1/13.