Modeling of Tumor-Immune System Interaction with Stochastic Hybrid Systems with Memory: A Piecewise Linear Approach

In this work, we bene t from hybrid systems that are advantageous because of their analytical and computational usefulness in the case of inferential modeling. In fact, many biological and physiological systems exhibit historical responses such that the system and its responses depend on the whole history rather than a combination of historical events. In this work, we use and improve hybrid systems with memory (HSM) in the subclass of piecewise linear di erential equations. We also include stochastic calculus to our model to exhibit uncertainties and random perturbations clearly, and we call this model stochastic hybrid systems with memory (SHSM). Finally, we choose tumor-immune system data from the literature and show that the model is capable to model history dependent behavior.


Introduction
Tumor dynamics exhibit complex interactions such as immune responses. Immune response to tumor growth has been widely investigated in the literature. Through these models, important parameters have been obtained and some predictions have been estimated and some predictions have been performed [30]. Some other approaches and models can be given as in references [5,7,12]. For a more detailed analysis of further models in the literature, one may refer to [1], and the references therein. In the case of potential use in treatment planning, instead of applying rst-principles models which require a lot of a priori information which might not be available, we employ a model that has the ability of adapting to a subject by subject variability and to unknown factors; this has the advantage of suggesting the best choice for each case. One of the well-known solution approaches to this problem is, inferring parameters and deciding on the behavior of the system based on empirical observations. In inferential modeling, a class of model is chosen and the parameters are derived from the observations. In such a case, hybrid systems, which can be described as systems formed by continuous and Boolean variables regulating each other [24], are useful in inferential modeling because of their analytical and computational advantages. Hybrid systems also oer several advances for various modeling [14,16] and theoretical problems [15] in natural sciences. In order to investigate dierent applications and various modeling versions and analyses of hybrid systems, one may see [6,19,32] and the references given there. Moreover, for piecewise linear approaches applied in regulatory systems, [37] oers a good source.
One of the most important properties of a regulatory system is that it is able to memorize its history. In other words, a combination of the previous inputs of the system decides its stationary behavior, and in turn, the system's stationary state decides on the response of the system's future external input. This is the crucial mechanism for adoption and learning in these systems. Delbruck's suggestion mentions that multistationarity, which is the existence of more than one stationary steady state of a system, might be reproduced by epigenetic dierentiation in gene networks [10,42]. Subclasses of hybrid systems, like piecewise linear ones, which are typical examples of non-trivial multistationary systems, are successfully used in modeling gene networks [14]. Some other models of similar phenomena which are also abstractions of hybrid systems, such as Boolean Networks [26,42] and Boolean Delay Equations (BDE) [9,35] are given, too. History-dependent responses of perturbations, or future inputs, can be obtained through the co-existence of multistationarity with delays. In [41], it has been stated that time delays could aect the stimulus responses in gene networks and, therefore, the answer of a model to a stimulus is history dependent. Some of such behaviors were demonstrated with the help of BDE [36]. By using hybrid systems with memory phenomena, one can eciently model the representation of factors which may depend on the whole history, rather than on a combination of historical events. Systems which use existing past information and involve it for intervention or decision making are known to be history dependent: associative memory (neural system) [3] and immune response [23]. Therefore, to consider functional extensions of hybrid systems is a strong motivation. Furthermore, involving delays in the piecewise constant part of hybrid systems was suggested to handle the analytical and computational diculties of traditional delay dierential equations (DDEs) with arbitrary initial functions [25]. Dierent types of Functional Dierential Equations (FDEs), which are generally developed by considering naturally occurring delays in dynamical systems [4], can establish the history-dependent behavior observed in biological systems [8,38]. For a more general and abstract model of hybrid systems with memory, one may refer to [31]. In this work, we construct a subclass of hybrid systems with memory to use in inferential modeling. In order to realize this goal, we propose to use Stochastic Hybrid System with Memory (SHSM) which diers from Hybrid Systems with Memory (HSM) by their ows and employ this approach to tumor-immune interaction dynamics. To nd the parameter values of the model, we benet from SDE Toolbox in MATLAB. We nd parameter values for every state of the model. Here, we benet from data that can be found the literature [13], in order to obtain more realistic values. The experiment that the authors of the reference [13] have conducted and the data values that they have obtained can be summarized as follows: Essentially, there are two dierent groups of mice and they are examined under IL1-α eect which is a variable in immune system. The group with IL1-α in their system is capable of defending itself against tumor growth, and the other group fails in this process. Their defense mechanisms and cancer growth processes dier widely. This shows us that the memorization ability of the immune system is crucial for the organism. The biological background can be found in Section 3. In this work, we consider a piecewise linear appraoch and for a nonlinear version, one may see [22].

Mathematical Model and Illustrative Examples
A typical hybrid system H = (Q, Y, Init, f, Inv, E, G, R) is formally dened as a combination of discrete states, continuous variables, initial states, vector elds, invariant sets, edges, guard conditions and reset maps [24,39]. Some extensions of hybrid systems are hybrid I/O systems which consider external inputs to the system [33], hybrid control systems involving a controller to determine a set of inputs forcing the system to supposedly optimal states [18], stochastic hybrid systems involving random perturbations to the system [29].
In the subsequent representation, a memory set has been included for the purpose of investigation of the behavior's dependence on history. Furthermore, we give the execution of the model. In fact, the behavior of the system and its response to external inputs are determined by the whole history rather than the initial values in many classes of switching systems in nature and technology [20,41]. Particularly, modeling some biological systems that exhibits history memorization it is necessary to involve an element to keep such information. More precisely, we have the following denition for Hybrid Systems with Memory.
A Hybrid system with memory H is [20,22]; consisting of • a set of discrete states Q = {q 1 , . . . , q m } which are also called locations, • a space of continuous variables X = R n , • a space of independent variables T = R k , typically the time T = [t 0 , ∞), • a vector eld f : Q × X × U × M −→ X, governing the continuous evolution, • an invariant set (domain, subspace) for each q ∈ Q, Inv : Q −→ P (X) where P (·) denotes the power set. Each state's governing dynamics is valid within its invariant set; • a set of edges (state transitions) E ⊂ Q × Q, • guard conditions for each edge G : • a reset map for each edge R : • M (t) ∈ M is a growing memory of past state transitions such that In this denition, the previous progress of the system, containing the values of variables and the time before and after the state transition, is sampled at state transitions. The memory set, M (t), is a piecewise constant variable between state transitions and, moreover, the memory grows at each state transition. Thus, a HSM has a complexity that is increasing with time. A subclass of HSM can be considered as a piecewise linear hybrid system with memory, and it can be described with a state space description as follows [20]: An example for the stochastic case can be given.
The system evolution at the initial state and until it hits one of the boundaries of the initial set b 1 or b 2 , looks as follows: The above stochastic system is Geometric Brownian Motion without drift; hence its solution can be found widely in the literature. We give the solution under Itô interpretation: The hitting times of each barrier are dened by stopping times as seen in the literature: Here, τ 1 and τ 2 are random variables. Thus, in each run, and when the rst transition occurs by hitting b 1 or b 2 is totally random. Thus, as a rst transition time we take τ * = min{τ 1 , τ 2 } and so the memory set will be Just as the system hits one of the boundaries, it experiences a dierentiation depending on the barrier that got hit.
Thus, our system can exhibit two dierent behaviors and, moreover, it has two dierent asymptotic distributions depending on its memory set. The stochastic systems used during each state are known as Ornstein-Uhlenbeck processes and have been widely used because of their mean reversion properties. Moreover, they have a stationary probability distribution unlike Brownian Motions. Thus, our system can develop according to two dierent stationary probability distributions and revert to two dierent means. The model solution is given by where y 0 is assumed as constant from the initial set. Thus, the distribution of the system after dierentiation is apparently a Gaussian distribution with possible two distinct means and variance values. The system's behavior will be where i = 1 if b 1 is achieved rst, and i = 2 otherwise.
3. Biological Background and Application on Tumor-Immune System Dynamics IL-1α inhibitors are being produced to avoid inammation, and the advancement of fever and sepsis; by this way, they treat diseases. In the regulation of the immune responses, they play an important role. For detailed information on IL-1α, one may refer to [11,13,34] and closer references provided there.
In [13], this eect of IL-1α inhibitors has been examined through some experiments, applied on dierent groups of mice. In a comparative study where the authors of [13] have investigated spleen cells from mice injected with an IL1-α-positive (Clone 2) cell line or a non-IL1-α-expressing (Clone 5) brosarcoma cell line, they have analyzed tumor-immune responses. All IL1-α-positive brosarcoma clones induced regression of tumor growth when they have been injected into mice. At rst stages of the experiment, cells started to multiply and accordingly regressed (see Figure 1). On the other hand, cells of non-IL1-α-expressing clones grow in a progressive manner and this results in the death of mice [13]. Both classes of cell lines (IL1-αpositive and -negative) exhibit almost identical growth progress, initially. Their work expresses the role of tumor cell-associated IL1-α in the induction of particular immune responses, ultimately leading to tumor regression and the improvement of an immune memory, which saves the mice from a battle with the violent tumor cells [13].
In our modeling framework, we will use the data of the experiment conducted by Dvorkin et al. [13]. We involve those data values which show tumor diameters of Clone 2 and Clone 5 according to days; in the setting of experiment Clone 2 has been injected with IL1-α, whereas Clone 5 has not been. As we have mentioned formerly, dierent levels of tumor growth and eector cells have been observed, depending on dierent values of IL1-α. Figure 1, Figure 2 and Figure 3 represent these eects graphically. Furthermore, the exact values may be observed from the Table 1. There, S.I. value refers to the Stimulator Index which corresponds to the ratio for immune cells (the eector cell and stimulator cells) and the tumor size has been measured in millimeters (mm) in that table. By examining the data values on the table, one can see that Clone 2 and Clone 5 behave similarly until day 3 and after day 3, Stimulation Index decreases in Clone 5, and after day 15 the tumor size is increasing in Clone 5. In order to illustrate those values in our model, we have designed a system in such a way that we partition it into 4 main states: rst state, i.e., q 1 , shows the behavior of both Clone 2 and Clone 5 until day 3; second state, i.e., q m , m = 11, 12, include slightly dierent behaviors of Clone 2 and Clone 5 until day 15; third and fourth states, i.e., q 2 and q 3 , respectively, show entirely dierent actions on tumor growth and S.I. values.
For SHSM simulation purposes, we have prefered Ornstein-Uhlenbeck type of stochastic dierential equations and splitted up the system due to the data values with the following manner: Depending on various memory values and dierent states, we have dierent parameter values: where q k ∈ {q 1 , q 11 , q 12 , q 2 , q 3 } and i = 1, 2, ... Network representation of the system can be seen in Figure 4. Now, let us consider a linear SDE: where α(t, X t ) = a 1 (t)X t + a 2 (t), For a general piecewise linear model, we use a similar statement in the works of Gebert et al. [17] and Öktem [37], with a version which includes the memory variables; and M s(t) 2 are matrices and k 1 and k 2 are vectors. Then the linear SDE can be represented by In this formulation, parameter values are thresholds and focal points. To determine these parameter values, one may refer to [37,40]. By exchanging equations into a two dimensional Ornstein-Uhlenberg process we get: or we can write in compact notation: In the equation (11), the term N q k ,m i is a focal point vector and M q k ,m i is a matrix. The memory set contains the rst transition in the rst state, time and state conditions. It will include and act according to q 11 or q 12 : where i = 1, 2. The described method below, that explains how to calculate the expected time for a stochastic process to reach a state or if a stochastic process can reach a state or not, can be found in [2,Chapter8]. By using the same process, we can calculate the time or the state of SHSM. In our model, we can nd those values compartmentally. For a stochastic process to reach states A or B, can be given by the following procedure. Suppose that X(t) be a stochastic process, and a solution of the stochastic dierential equation given by then the transition probability density function for the linear stochastic process gives us a solution to the backward Kolmogorov dierential equation [2,Chapter 8]: If Q(x, t) is the probability that measures whether the process does not reach states A or B in time [0, t], and if A < x < B, then we have [2,Chapter 8]: If T (x) is the random variable which refers the time for the stochastic process to reach states A or B, then expected time can be found by the following integral [2,Chapter 8]: In order to simulate the stochastic dynamics, and estimate parameter values, we have beneted from SDE toolbox in MATLAB. Besides the parameter values, we are able to nd Monte-Carlo statistics results including process mean, process variance, process median, condence interval for the trajectories, process skewness, process kurtosis, process moments by utilizing this toolbox (cf. Figures 5-9). Since we can nd the values of the unknown parameters, the equations for each state can also be determined. Some data values have been kept and modeled according to these data in order to test the model. You may see the Figures 10 and 11.      In that state, we expect the system to show a similar behavior as Clone 5. Variable X1 (left) and Variable X2 (right) [21]. In that state, we expect the system to show a similar behavior as Clone 2. The variable X1 (left) and the variable X2 (right) [21].

Discussion
The biological example used here can be thought of, e.g., as an example to an alternative way to use FDEs or DDEs. In the rst state, Clone 2 and Clone 5 acts similarly, and after day 3, the tumor size starts to change slightly in dierent clones. Clone 2 remembers the values of t, q 1 , x and acts between days 3-15 as it did in the past. Therefore, q m , which includes q 11 and q 12 and the host choosing one of the states according to its memory set, can be considered as the delay state for Clone 2 in the experiment. After day 15 it changes its behavior. This behavior can be considered in such a way that the immune system remembers the memory set of t, q, x and after day 15 it changes its behavior according to IL1-α which can be considered as the control input u of the host. The simulations have been done compartmentally.

Conclusions
In this work, we utilized Stochastic Hybrid Systems with Memory (SHSM) to model the history dependent behavior in hybrid systems and stochastic hybrid systems. The work can be considered either as introducing memory into stochastic hybrid systems or introducing switching behavior into FDEs. The introduced class can be used in modeling of many systems with pattern memorization capability and subject to external inputs. Various hybrid systems involving delays have already been investigated. We selected a system exhibiting the same pattern, which is simple but not easy to model with traditional hybrid systems. Various types of functional dierential equations were also implied in modeling history dependent systems by other scientists. We discussed stochastic case, where stochastic hybrid systems with memory can be a tool for modeling non-Markovian Stochastic Dynamical Systems. In this work, we have seen that our modeling scheme can be employed to model tumor-immune dynamics. We have worked with raw data found from the literature and have seen that complex networks, which exhibits history-dependent behavior, can be modeled in a simpler way by applying this modeling framework where the dynamics of the system is determined by the location of the state vector and the memory. For instance, gene-regulatory networks that has memorization capability can be mimicked by this approach. As a future work, SHSM can be employed in modeling some complicated genetic or neurological systems. Another possible use can be the design of learning adaptive systems or controllers. Moreover, since fractional dierential equations are said to capture and exhibit memory dependent behavior, they are used in modeling real world applications [27] and they are used to check stability where the equation is very sensitive to delays [28]. Therefore, an extension of this method to fractional dierential equations and checking and comparing the stability of the system can also be considered as a future work.