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, mechanistics 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] 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, resulting 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.

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 ( 2𝜋𝑅 𝑇 𝑠𝑖𝑛 8𝜋𝑡 3𝑇 − 𝐸𝑐𝑜𝑠 8𝜋𝑡 3𝑇 ) + 𝐸𝐴 2 + 𝐸𝐵 if 𝑡 ∈ (0, 3𝑇 8 ] 𝐴 2 ( 2𝜋𝑅 𝑇 𝑠𝑖𝑛( 8𝑡 5𝑇 + 2 5 )𝜋 − 𝐸𝑐𝑜𝑠( 8𝑡 5𝑇 + 2 5 )𝜋) + 𝐸𝐴 2 + 𝐸𝐵 if 𝑡 ∈ ( 3𝑇 8 , 𝑇 ](1)

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𝑇 8 ] −𝐷𝑜(𝑡) if 𝑡 ∈ ( 3𝑇 8 , 𝑇 ](4)

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.

PID Controller & Simulations

To begin with, a feedback controller will be added to study different respiratory control mechanisms. Similarly to [10], who inlude a closed-loop respiratory system model of the brainstem 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 development 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. 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 each instant, there is a target 𝑜 * 𝑎𝑣 , and the controller attempts to achieve it by modifying either 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 target value 𝑜 * 𝑎𝑣 . Table 1 includes 9 different changes in 𝑜 * 𝑎𝑣 , one in each row. A change occurs 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 𝑜 * 𝑎𝑣 .

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 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𝑘.

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.

𝑢 𝑓 𝑓 𝑡 = 𝑔(𝑥 * 𝑡+1 , 𝑥 𝑡 ; w)(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.

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.

Immune equations and parameters

Cytokine storm (section 5.1) is a life-threatening inflammatory response characterized by hyperactivation 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

𝑑𝑉 𝑑𝑡 = 𝑣 𝑖𝑛 𝐼 (𝑡) − 𝑘 𝑣 𝑉 (𝑡)(10)

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 +𝑥(𝑡) (𝑥(𝑡) − 𝑚 𝑖𝑠 )(𝑦 1 − 𝑦(𝑡))(𝑦(𝑡) − 𝑦 2 ). The normal transformation of cytokines is given by 𝑦 𝑖𝑛 − 𝑘 2𝑖𝑠 𝑦(𝑡) and the production stimulated by immune cells 𝑏 2 𝑎 1 𝑦(𝑡)𝑥(𝑡) 𝑎 2 (𝑐 2 +𝑥(𝑡)) . Thus, the followings equations have been reproduced in Figures 3 :

𝑑𝑥 𝑑𝑡 = 𝑥 𝑖𝑛 − 𝑘 1𝑖𝑠 𝑥(𝑡) + 𝛾 1𝑖𝑠 𝑥(𝑡)𝐼 (𝑡) + 𝑏 1 𝑥(𝑡) 𝑐 1 + 𝑥(𝑡) (𝑥(𝑡) − 𝑚 𝑖𝑠 )(𝑦 1 − 𝑦(𝑡))(𝑦(𝑡) − 𝑦 2 ))(12)

𝑑𝑦 𝑑𝑡 = 𝑦 𝑖𝑛 − 𝑘 2𝑖𝑠 𝑦(𝑡) + 𝑏 2 𝑎 1 𝑦(𝑡)𝑥(𝑡) 𝑎 2 (𝑐 2 + 𝑥(𝑡))

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 ) = 𝐷(𝑡 𝑖 ) + 0.1 𝑦(𝑡 𝑖+1 ) − 𝑦(𝑡 𝑖 ) 3(14)

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. 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.

𝐼 (𝑡𝑖+1 ) = 𝐼 (𝑡 𝑖 ) − 0.1(𝑦(𝑡 𝑖+1 ) − 𝑦(𝑡 𝑖 ))(15)

Conclusions

In summary, the variables that have been simulated for the respiratory system are: CPG amplitude 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.