Simulations⋆

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 incorporating 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 diferent medical systems. Thus, helping to bridge the gap towards more explainable 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] EXPLIMED - First Workshop on Explainable Artificial Intelligence for the medical domain - 19-20 October 2024, Santiago ∗Corresponding author.

-

CEUR ceur-ws.org

1. Introduction

emphasize, causality is indeed crucial for explanatory purposes since, thanks to it, an efect 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 diferent 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 diferential 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 efect 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 efective 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 ( 2 ( 2 2 ( 8 3

− 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 ()

in the instant of time is: the amplitude. Furthermore, it is considered that the respiration is divided in two phases: inspiration occurs when ∈ (0, 38 ) and expiration occurs when ≥ 3 . The muscular activity 8 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() +

2() ′() = 1 [− () + ()] Finally, the oxygen concentration () at time is: ′() = { −() + ( −() − ) ′() if ∈ (0, 3 ] if ∈ ( 3 , ]

8 8 where = 0.75 −1 is the difusion 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 diferent respiratory control mechanisms. Similarly to [ 10 ], who inlude a closed-loop respiratory system model of the brainstem (1) (2) (3) (4) 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 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 diference between the current value for the variable versus the desired value for the same variable. The PID control algorithm consists of three diferent 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 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 diferent changes in ∗ , one in each row. A change occurs every 10 seconds as displayed in Figure 2, and each change has a diferent 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 ∗ .

4. Feedforward Control & Inverse model

A PID controller is a control mechanism that, through a feedback loop, allows to calculate the diference 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 .

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: (5) (6) +1 = ( , ) = ( , + ) 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)

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 afects the amplitude , or period , or both. This, in turn, changes the CPG control signal () on the muscle () , which also afects 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

section we will describe a linear inverse model: 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 (7) (8)

= 1 ⋅ + 2 ⋅ ∗ 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, ∇ =

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 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 diferential equations. [ 16] applies a simple two-component diferential 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 fiteen ordinary diferential 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:

= − () − ()() = ()() −

() − () () = () − () (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+( )() (() − )( 1 − ())( () − 2). The normal transformation of cytokines is given by − 2 () and the production stimulated by immune cells 2 21((2)+( )()) equations have been reproduced in Figures 3 : . Thus, the followings = − 1 () + 1 () () + (() − )( 1 − ())( () − 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 efect 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 difusion rate of oxygen concentration and the rate of increasing oxygen concentration , that afect the amplitude, period, or both, to achieve the target .

() 1 1 + () = − 2 () + 2

1 ()() 2( 2 + ()) ( +1 ) = ( ) + 0.1 ( +1 ) − ( ) 3

(9) (10) (12) (13) (14)

( +1 ) = ( ) − 0.1( ( +1 ) − ( )) (15)

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 dificult, 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 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 efect 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 efective therapeutic strategies. Thus, this is a first step in simulating the efect of diferent pathologies, and it will be the basis to also simulate the efect of diferent (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. [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.

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 linear 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 focused 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 dynamics during a cytokine storm , Mathematical and Computational Approaches in Advancing Modern Science and Engineering ( 2016 ) 331 - 339 . doi: 10 .1007/978- 3- 319 - 30379- 6_ 31 .