Stochastic approach to study the properties of the complex patterns observed in cytokine and T-cells interaction process

Patterns in complex systems store hidden information of the system which is needed to be explored. We present a simple model of cytokine and T-cells interaction and studied the model within stochastic framework by constructing Master equation of the system and solving it. The solved probability distribution function of the model show classical Poisson pattern in the large population limit $M,Z\rightarrow large$ indicating the system has the tendency to attract a large number small-scale random processes of the cytokine population towards the basin of attraction of the system by segregating from nonrandom processes. Further, in the large $\langle Z\rangle$ limit, the pattern transform to classical Normal pattern, where, uncorrelated small-scale fluctuations are wiped out to form a regular but memoryless spatiotemporal aggregated pattern. The estimated noise using Fano factor shows clearly that the cytokine dynamics is noise induced process driving the system far away from equilibrium.


I. INTRODUCTION
Patterns are inherent fundamental building blocks of complex systems.Each pattern is created by a self-organized certain population of interacting components of the system resulting aggregation having certain form of regularity, symmetry and structure which is the source of certain form of information and can participate in the system's collective decision making process [1].However, the whole process of each pattern may not be the exact sum of the individual activities of the constituting components which are nonlinear in characters, sensitive to the surrounding internal and external fluctuations and initial conditions, dissipative in nature, and the dynamics are generally strongly coupled [2].Each pattern is generally open in nature which allow exchange of energy, accessible states are far away from equilibrium, and the information flow is permitted across the boundary [3].The interaction among the components in a pattern could be of diverse form which drives the dynamics of individual components quite complicated in nature which could be the origin of complicated properties of each pattern [4].Technically, the degree of pattern is generally characterized by complexity [5].The complexity of the pattern can be measured by different techniques, some of which are information theoretical approach [7], statistical complexity measures [6], Kolmogorov complexity [20], Lempel-Ziv complexity [9] etc.The study of the properties of these patterns may provide the information of how complex systems work at fundamental level.
The observation of the emergent properties of the patterns in the systems can be done mathematically using traditional central limit theorem by calculating large number of outcomes perceived from hidden processes of the distribution of the interacting random variables {x i }; i = 1, 2, ..., N corresponding to the aggregated N components of the pattern [11][12][13].Depending on the observed properties of the probability distribution, the pattern can be in one of the classical patterns, namely, Norma, Poisson, Power law patterns etc [10] and these patterns, in most of the biological systems, are dynamic in nature depending on the fluctuations driving them and limiting distribution conditions [14].If the spatiotemporal events are scattered randomly throughout the system, the probability distribution function becomes Poisson distribution and the pattern corresponds to classical Poisson pattern [10,12,13].Since Poisson pattern is generally generated from the random processes smaller-scales, the basin of attraction of such patterns have aggregation of large scale small-scale random processes by segregating from nonrandom processes [15].On the other hand, classical Normal pattern stores the information of mean x and variance σ 2 in the aggregation process of the components in the pattern P[ x , σ 2 ], where, sum of the uncorrelated small-scaled fluctuations of the components during the aggregation process cancels out to attract towards a certain form of spatiotemporal regular aggregated pattern [10,11].Hence, the dynamics of the aggregated components in such classical Normal patterns follow Brownian motion [18].Further, if the aggregated pattern keeps only the information of mean, the pattern falls in exponential pattern P , where, the pattern is contributed from the waiting time characterization of random events and the distribution has memoryless property [10].If the aggregation pattern preserves the information of geometric mean, the pattern follows classical Power law pattern which can have fractal patterns which have self-similar process , where, h is self-similarity exponent [16,17].From these probability distribution functions, one can measure hidden information stored in the respective patterns by calculating Shannon entropy, The time evolution of the components in the patterns/systems is in general fluctuation driven stochastic process which can be well described in a probabilistic manner [20].The inherent fluctuations in the dynamics of the constituting variables arise from the internal (due to inter-molecular interaction) and external (due to fluctuations driven by surrounding the system i.e. temperature fluctuation, environmental fluctuations etc) fluctuations and have variety of roles in regulating the system [21,22].The trajectory of the variables of the pattern/system and their realizations can be obtained using the Master equation, which is the discrete form of the Chapman-Kolmogorov equation [20,26], constructed from the information of the transitions of the states due to inter-variables interaction [22][23][24].On the other hand, computationally, the trajectory of each variables in the system can be traced using Monte carlo like stochastic simulation algorithm (SSA) due to Gillespie [25].The stochastic approach has wide variety of applications, ranging from microscopic biomolecules [27] to macroscopic ecology [28] to cosmology [29], and then to finance [30,31], cryptography [32] etc.
Cytokines are certain category of diverse family of small proteins or glycoproteins that are secreted by specific cells of immune system and are produced throughout the body by cells of diverse embryological origin [33].Cytokines provide the mechanism by which lymphocytes, inflammatory cells and haematopoietic cells can communicate with each other to mount and coordinate an effective immune response [34].Interleukins (ILs) are the Cytokines secreted by leukocytes and affect the cellular responses of leucocytes [55].Cytokines are pleiotropic which means that they have multiple effects and can increase both proliferation and death of a certain cell type [36].On the other hand, a pleiotropic cytokine that helps in the growth of T-cells is Interleukin-2 (IL-2) and intervenes in the activation-induced cell death (AICD) and causes differentiation of regulatory T-cells [37].Hence, interleukin-2 (IL-2) plays a major role in the proliferation of activated T cells [38].These T-cells can identify and eliminate invading pathogens, and also attack human tissues [39].Further, a successful activation and differentiation of T-cells is aided by the presence of cytokines and specify the type of response that the T-cells give to a foreign antigen [40].Now, it has been shown that, a particular type of T-cells known as the CD4 + T-cell produces cytokine IL-2, which is responsible for both the proliferation and death of T-cells [41,42].This cytokine can also be produced by T-cells, the T-cells proliferate to the cytokine again and are removed, they can also uptake and secrete the cytokine and the cytokine can then degrade itself [36].However, in the complicated regulating mechanism of cytokine and T-cells interaction what are dynamical properties of the system, how does the dynamics affect the patterns in the interacting system and how can we interrelate the emergent patterns to phenotyical behavior are still not answered fully.Further, what are the roles of the fluctuations in system's dynamics are still open questions.
In this work, we present a simple model of cytokine and T-cells interaction model.We then analyze the model using stochastic approach to understand how the probability distribution function of the cytokines aggregate to form different patterns at various situations.The results also explain the role of noise in regulating the system.Then some conclusions are made based on the results obtained.

II. METHODOLOGY A. Master Equation Formalism
The dynamics of stochastic particles can generally be described by The Chapman-Kolmogorov equation which is based on the probabilistic description of state transitions during the course of the particle's motion [20,26].This Chapman-Kolmogorov equation in integral form [26] The standard way of deriving Chapman-Kolmogorov equation in differential form is to calculate the time evolution of the expectation value of a function f (x) defined in the n−dimensional system given by, Now, using the definition of differentiation to the partial derivative of the probability distribution function P (x, t|x ′ , t ′ ), expanding f (x) with Taylor series expansion and then rearranging the terms by keeping terms upto second order derivative, one can get the following Chapman-Kolmogorov equation [20,26], where, {W } are Wiener functions which are transition probabilities of the states and A i and B ij are the functions which depend the properties of the system.The continuum approximation of the equation ( 3), which is contributed from the first two terms in the right hand side, provides Fokker-Planck-Kolmogorov equation [26].The Master equation is generally the discrete form of the Chapman-Kolmogorov equation [20,22,26] which is contributed from the jump process of the particle dynamics (keeping the last two terms in the right hand side), which is given by, This equation can also be expressed by replacing integration by summation over configurational states indicating jump process [22,26,27] which is easier to handle and more directly related to physical systems [22].Master equation can also be seen as the time evolution of the configurational probability for Markov processes of systems that jump from one to another state which is due to birth and death processes of the particles in the system.There are different techniques to solve Master equation, namely, using generating function technique [24], Fourier transform technique [29], Laplace transform technique [48], Laplace-Fourier transform technique [49] etc.However, solving Master equation of complex systems involving multi-variable is quite difficult except for simple systems.

B. Generating function technique to solve Master equation
The generating function (GF) technique is generally used to solve Master equations for simple low dimensional systems [24].Brief procedure to solve Master equation using GF is given as follows.The n−dimensional GF G(s 1 , s 2 , ..., s n ) is the transformation map of the probability distribution function P (x 1 , x 2 , ..., x n ; t) as defined by, where, {N [0] i }; i = 1, 2, ..., n is the set of initial populations of the respective variables x 1 , x 2 , ...x n .Now, the solution G is put back to the GF (5) and equating the coefficients of the terms containing n i=1 s xi i from both sides, one can arrive at the functional form of P (x 1 , x 2 , ..., x n ) which is the solution of Master equation (4).However, solving Master equation for any high dimensional system is quite difficult and is still an open question.
The generating function also allows us to determine the observables which will describe the properties of the system, such as, the role of fluctuations in the system and many others.The analysis of GF and probability distribution function could reveal various other hidden patterns and their properties of the system.
The observables of the system namely, mean, variance and other measurable parameters can be determined using the generating function.The mean of any variable x i can be obtained from the following relation, x The fluctuations involved in the dynamics of the particles in the system can be well quantified by measuring Fano factor F of the variable corresponding to the particle which is given by [50], The Fano Factor is an important parameter which can determine the role of fluctuations or noise in the dynamics of the particles of the system [51].Depending on the values of F , one can determine the behavior of the particle which undergo various processes [51,52].
• If F = 1, then, the process is a Poisson process and the particles in the system undergo Brownian motion.In such situation, there is no specific role of the fluctuations in the dynamics of the system.
• If F > 1, the process is a noise-enhanced process and is also known as super-Poisson process.Here, the fluctuations has significant role in regulating the dynamics of the system.
• If F < 1, the process is called a sub-Poisson process.Here, the noise has no significance but still has a tendency to regulate the dynamics of the system.

C. Cytokines and T-cell interaction model
The proliferation of T-cells (T) involves the production of a cytokine IL − 2 (Z), and this cytokine plays a significant role in the T-cell proliferation and death [36].In this simple model, if we denote the rate of proliferation as X and rate of death as Y , and when these two rates are functions of the produced cytokine Z, we can write the concentration of the cytokine Z at any instant of time by a simple first order differential equation [36], The maximum value of Z during the process can be obtained by calculating the fixed point of the equation (11), T .This steady state population of the cytokine Z is generally independent of its initial population at time t = 0.However, the initial population of cytokine trigger the T-cell proliferation by producing another type of cytokine in response [36].Moreover, the cytokine IL-2 can be produced from other external sources other than T-cell proliferation.Further, T-cells can uptake and secrete the cytokine IL-2 which are needed to be incorporated to the model.
Considering various processes involved in the proliferation of T-cells and its death with the production of cytokine IL-2, the above model (11) can be represented by following form of the equation [36], where k 1 is the secretion rate by an external source, k 2 is the secretion rate by T-cells, k 3 is the rate of uptake of Z by T-cells and k 4 is the rate of degradation of Z by itself.This equation can translated into the following four reactions,

A. Stochastic approach to study cytokine and T-cells interaction model
The study of the stochastic dynamics of the system can lead us to explore the role of fluctuations and other factors which regulate the behavior of the system.Since the patterns exhibited in the interaction of the cytokine and T-cells at any instant of time is random and probabilistic in nature, understanding the properties of the patterns within the stochastic framework could be important to correlate with the exhibited phenotypic behavior.We consider the equivalent reaction set (13) of the model equation (12) and assume that the total population of cytokine Z and T-cells is a constant, i.e, T + Z = constant = M .Further, as the T-cells recruits cytokine in a quite fast rate [53] and the T cells population in our immune system always tends to an equilibrium/steady state [54], the population of T-cells is also taken to be a constant as compared to cytokine population in our calculation.
Consider the transition of states given by the model reaction set (12) obey the Markov property [33], which means that the future state of a system depends only on the present state and not on any of the past states, we can determine the transition probabilities from one state to the next for each reaction using the Law of mass action.From the model (12), the reaction α k1 −→ Z indicates the production of cytokine Z by T-cells, where, α is a constant source.Since the firing of the reaction follow Markov process, the transition rates of the process can be calculated using mass action law in the state change processes and are given by, ω 1 = αk 1 ∆t + O(∆t 2 ) and ω ′ 1 = αk 1 ∆t + O(∆t 2 ) respectively.Next, the reaction T k2 −→ Z reveals the proliferation of T-cells to recruit cytokine Z and removal of itself with a rate k 2 .The corresponding transition rates are given by, ) respectively, where, β is the mean value of population of T-cells as its fluctuation is approximated to be negligible.Further, the reaction T + Z k3 −→ T shows that the T-cells can also uptake a cytokine Z, secrete another cytokine Z and then remove itself.This means that the state of cytokine Z remains constant at state Z.The transition rates of this reaction are given by, ω 3 = k 3 β(Z + 1)∆t + O(∆t 2 and ω ′ 3 = k 3 βZ∆t + O(∆t 2 respectively.Lastly, the degradation of cytokine Z is indicated by the reaction Z k4 −→ φ.The transition rates are given by, ω 4 = k 4 (Z + 1)∆t + O(∆t 2 ) and ω ′ 4 = k 4 Z∆t + O(∆t 2 ) respectively.Now, according to the principle of detailed balance, the process of transition from one state to another should be in equilibrium with the corresponding reverse process of transition.Using this principle, we can determine the probability, P (T, Z; t + ∆t), that the system is in the state (T, Z) at any instant of time t + ∆t provided that the system was in some other state at time t.Now, using the transition rates for each reaction, we can arrive at the following detailed balance condition, Now, rearranging the above equation and using the definition of differentiation in the limit ∆t → 0, we arrived at the following Master equation, This Master equation ( 15) can be solved using generating function technique which we have explained in the methodology section.The generating function in one dimension can be defined as.
with the initial condition, G(s, 0) = S N at t = 0, where N is the original number of cytokine population in time t = 0. Now, multiply equation ( 15) by s Z on both sides, summing over Z, and after rearranging the terms, we have, This spatiotemporal partial differential equation ( 17) can be solved using standard Lagrange's Charpit characteristic method for partial differential equations [55].Now, applying this method, the solution can be obtained from the following set of equations, Now, let us consider the first equation constructed from the above set of equations ( 18), dt = ds (βk3+k4)(s−1) .Integrating both sides and rearranging the terms, we get, where, λ 1 is a constant and C 1 is the constant of integration.Next, we consider second equation from the set of the equations eqreflce, given by, . Now integrating and rearranging the terms, we have, where, λ 2 = e C2 is a constant and C 2 is another constant of integration.Now, combining the equations ( 19) and ( 20 Now, to find out the functional form of the function F in equation ( 21), we take the initial boundary condition given by, where, M is the initial cytokine population at t = 0. Putting this condition to equation ( 21), we have the F given by, Then, using this equation to the equation ( 21), we have the GF given by, where, a = e −(βk3+k4)t and b = αk1+βk2 βk3+k4 .Now, rearranging the first factor of the equation ( 24) in the form F 1 = (s − 1)e −(βk3+k4)t + 1 M = a M s + 1−a a M , and using Binomial expansion, we have, The second factor in the equation ( 24) can be rearranged and expand in the form, putting the factors back to the equation ( 24), and taking Z = i + j, we have, Equating the coefficient of s Z , we arrived at the expression for P (Z, t), given by, Now, let us study the role of noise in cytokine dynamics triggered by T-cells which can be analyzed by calculating Fano factor [50].To calculate the Fano factor F of the cytokine dynamics we need to calculate observables, mean Z and variance σ 2 = Z 2 − Z 2 .Using equation eqrefgf, the mean of the cytokine dynamics is given by, The asymptotic behaviors of the mean can be calculated as follows.Taking lim t→0 Z = e b M , which indicates that the initial population of cytokine Z depends on various interaction processes involved in the model i.e. various rate constants and population of T-cells β.Since the total population has to be M , we should have b = 0, which means that αk1+βk2 βk3+k4 = 0 i.e. αk 1 + βk 2 = 0 revealing that the overall production of cytokine both from external sources and secretion by T-cell should be zero.On the other hand, since, lim t→∞ Z → be b = constant indicating the population of cytokine reach a stable constant value at significantly large time.Next, the variance is given by, Now, the Fano factor is calculated as follows, Now, we analyze the Fano factor to understand the behavior of the cytokine dynamics.

Proposition 1
The asymptotic behavior of the cytokine dynamics allows only sub-Poissonian process.
Proof: Let us calculate the asymptotic values of Fano factor F from equation (29).
From the expression of b, b = 0 and b 0, such that, e b 1.Hence, lim t→0 F → −ve which can not be physically reliable and can not say anything about the cytokine dynamics.On the other hand, This indicates that the cytokine process follows sub-Poissonian nature.The involved noise in the cytokine dynamics tries to stabilize its dynamical behavior.
Proposition 2 At any instant of time t 0 t ∞, the cytokine dynamics noise-enhanced process.
Proof: The cytokine dynamics becomes noise enhanced process or super-Poissonian process when F satisfies F > 1.
Applying this condition in equation ( 29), the condition becomes, Since, 1 ≤ e b < ∞, 1 − e b < 0. This indicates that left hand side term of equation ( 32) is negative value.Since, Z 2 > 0, the condition (32) is always satisfied.Further, since Z 2 can not be negative the process can not be Poissonian and sub-Poissonian process.Hence, only possible role of noise for 0 t ∞ is noise enhanced process.

B. Patterns observed in cytokines distribution dynamics
The simple cytokine and T-cells interaction model exhibits complicated system's dynamics where the role of noise is quite significant in regulating the dynamics.Interesting question is even though the dynamics of the cytokine molecules are complicated they should follow some traditional equation of motion.Since the typical size of the single cytokine molecule (IL2) is around 450nm in diameter [56] which is within mesoscopic scale size [57], they are very sensitive and significantly affected by internal and external fluctuations [58].We intend to study the patterns of the distribution functions these particles by analyzing the probability distribution function given by equation ( 26).
Theorem 1 Near thermodynamic limit, M, Z → ∞ and Z → f inite, • The cytokine distribution function P (Z, t) follows modified Poisson distribution: which is independent of cytokine population variable Z.
• The information stored in classical Poisson pattern is given by, Proof: At large M limit, the equation ( 26) can be written as, The Poisson distribution is independent of Z indicating the cytokines' distribution becomes homogeneous approximately everywhere in the system when the population of Z is significantly large.Now, the complexity in the distribution of the cytokine molecules whose dynamics undergo stochastic process can be calculated using Shannon entropy H using the modified Poisson probability distribution function P(Γ, t) = P ois(Γ, t) given by equation (35), where, Γ = be a( M b −1) .We follow the procedure of calculation given in the work [43].Then the terms are rearranged to get the expression for H as in the following, Then substituting the series expression for ln(i!) given by the Malmsten's representation [44] as in the following, The equation (36) becomes, The integration in equation (38) can be done using series of formulae listed [44], doing integration by parts and rearranging the terms, we get the following expression for entropy [43], Now, let us calculate the information stored in the classical normal pattern which can be calculated using the equation (40).The Shannon entropy, which can be used to calculate the complexity in complex system, can be obtained by using the equation (43) and rearranging the terms as given below, We, then, substitute the expression for Normal distribution given by N .Then, after integration and rearranging the terms, we have, where,

√
π are constants respectively.erf (x) is the error function of x.From equations (40) and (46) and putting the value of H max = e −a(M+b)+Γ × ln(M ), we have the information stored in classical Normal pattern, I as in the following, In order to have non-negative information, one has to have the condition ln . We further get that at the sufficiently large time, . This indicates that at large time limit, the overall production rate of cytokine should be controlled such that ln(M ) ≫ b as the third term is less than one, otherwise, I ∞ < 0 which is physically have no meaning.

IV. CONCLUSION AND DISCUSSION
Patterns are manifestation of regular self-organised aggregated random components.They have significant roles as they keep diverse information stored in them to take part in collective decision making.Hence, it is important to study these patterns and information associated with them specially in biological systems.
The interaction of cytokine and T-cells in our body keeps the regulation of immune system and helps in various biological functions [46].Since the cytokine production by T-cells is quite essential for developing immune system, we considered a simple of cytokines and T cells interaction, where, all possible cytokine production mechanisms are associated.We used stochastic formalism to study the model by constructing Master equation of the model system specifically considering the dynamics of cytokines.The constructed Master equation is solved analytically using generating function technique and the solved probability distribution function, which describes the pattern of distribution of dynamical cytokines in the system, is found to be quite complex.However, under certain limiting conditions, we explored few interesting patterns which carry the information of the system.Near thermodynamics limit, the probability distribution function becomes modified Poisson distribution which indicates that the distribution of the cytokines in the system exhibit classical Poisson pattern.This property indicates that the cytokine events occurring in the system is purely random and the pattern is exhibited mostly by the self-organised aggregation of spatiotemporal random cytokine events scattered over the entire system.In such situation, Poisson pattern behaves as a basin of attraction which attracts cytokine random processes towards it and nonrandom processes are away from it [10].The uncertainty or surprise of random cytokine events is measured by complexity measurement parameter Shannon entropy which we found to be time dependent.The information stored in the Poisson pattern is then calculated as the difference of maximum entropy and entropy at any instant of time, and found that the competition between overall production rate of cytokines and overall decay rate of it need to be balanced such that information stored in the Poisson pattern should be optimised and non-zero.
On the other hand, near the thermodynamics limit, if the mean of the cytokine distribution is taken to be significantly large, then the distribution of the cytokines exhibits classical Normal pattern.In such situation, the phenomenon indicates that the sum of the observations of the small-scale random and independent cytokine events is also a random variable which is a larger-scale random process in the system, no matter how small the distributed random variables are [47].The distribution of such larger-scale random processes in the system approaches to classical Normal pattern.It also reveals that not only the small-scale random processes, but also larger-scale random processes of the cytokine distribution exhibit Normal pattern.This could be an indication of hierarchical organization of the cytokine random processes which approaches to similar behaviour at various scales.The information contained in the Normal pattern is calculated using Shannon entropy of the distribution.
The role of the noise is important and significant specially in microscopic systems.In the cytokine and T-cells interaction model, the strength of the noise is measured using Fano factor which can be calculated from the information of mean and variance of the distribution.The calculated Fano factor shows noise enhanced or super-Poissonian process which indicates that the noise associated with the distribution can significantly drive the system at various possible accessible states.In such situation, the state of the cytokine distribution and dynamics is driven far away from the equilibrium and are generally in active state.
The study in this simple model could able to understand few important features and properties of cytokines distribution which is dynamic in principle in cytokine and T-cells interaction.This system is multi-variable system which is needed to incorporate various other interacting molecular species.In such situation, analytical work could have limitation but definitely large scale simulation can be implemented to the system to study dynamics and distribution of the cytokine and related other variables.
) Now, multiplying the above Master equation (4) by x1,x2,...,xn n i=1 s xi i and substituting equation (5), one can transform the Master equation to spatiotemporal partial differential equation in the GF.The partial differential equation in G can be solved using the boundary condition at t = 0,

45 )
to obtain the expression for H.At the M → ∞ limit, can relate the joint probability distributions of different sets of variables involved in a stochastic process where the transition of states of the system obey Markov process.If one consider a ′ n ′ -step transition probability of the stochastic process, then the jump probability P (x n , t n |x 1 , t 1 ) from the state (x 1 , t 1 ) to the state (x n , t n ) is given by, P (x n , t n |x 1 , t 1 ) = ... dx n−1 ...dx 3 dx 2 P (x n , t n |x n−1 , t n−1 |x n−2 , t n−2 )...P (x 2 , t 2 |x 1 , t 1 )