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.