Controllability in hybrid kinetic equations modeling nonequilibrium multicellular systems.

This paper is concerned with the derivation of hybrid kinetic partial integrodifferential equations that can be proposed for the mathematical modeling of multicellular systems subjected to external force fields and characterized by nonconservative interactions. In order to prevent an uncontrolled time evolution of the moments of the solution, a control operator is introduced which is based on the Gaussian thermostat. Specifically, the analysis shows that the moments are solution of a Riccati-type differential equation.

1. Introduction

The derivation of mathematical models for the evolution of multicellular systems is a challenge of the last and also this century. Several mathematical frameworks have been proposed with the aim of capturing the main properties of these systems, which have also the feature to be complex [1]. Indeed different emerging behaviors appear as result of the interactions: competition, pattern formation, and organization (see the book [2] and the review paper [3]).

The mathematical frameworks present in the pertinent literature principally pertain to multicellular (biological) systems that are not subjected to external forces (systems at equilibrium); the interested reader is referred to the recent review [4] (and the references cited therein) for what concern the hyperbolic and kinetic models for living systems at equilibrium.

However, dealing with multicellular systems whose time evolution depends on the presence of external fields requires a different treatment. An external force field moves the system away from equilibrium and the blow-up of some observable (such as the kinetic energy) occurs.

Recently a new mathematical framework has been proposed and analyzed: the thermostated kinetic theory for active particles [5-9]. This framework, which refers to general complex systems (including vehicular traffic, and crowd and swarm dynamics; see, among others, [10-13]), is based on kinetic theory for active particles and deterministic thermostat. Kinetic equations coupled with the Gaussian isokinetic thermostat are already present in the literature: Kac equation [14-16] and Boltzmann equation with the one-dimensional BGK operator [17] (see also the recent review paper [18]). The solution of these equations is a distribution (probability) function defined in the microscopic state of the particles, which can include mechanical variables (space and velocity) and/or biological variables.

Roughly speaking, the deterministic thermostat is a mathematical framework that has been proposed in nonequilibrium molecular dynamics in order to achieve equilibration to nonequilibrium situation where there is a flux of energy through the system induced by external fields. The reader interested in further details is referred to papers [19, 20] and the review paper [21]. Accordingly, a control term (a damping term) is inserted into the equations of motion in order to ensure the conservation of the kinetic energy (or a different moment); see [22, 23].

The thermostated kinetic for active particles framework proposed in [7] is concerned with nonequilibrium complex systems characterized by conservative interactions only (namely, interactions that only modify the microscopic state of the particles). Therefore, the mathematical modeling of nonequilibrium multicellular systems, where proliferative and destructive events may occur, cannot directly performed within this framework.

This paper is concerned with a further generalization of the thermostated kinetic framework [7] that takes into account the modeling of multicellular systems with nonconservative interactions. Moreover, the velocity variable is discretized, generating an hybrid framework with discrete and continuous variables. Therefore, the new thermostated kinetic framework consists of nonlinear hybrid partial integrodifferential equations with quadratic nonlinearity.

It is worth stressing that, to the best of our knowledge, this is the first time that nonconservative interactions are taken into account in the thermostated kinetic for active particles framework.

Specifically, in the present paper the time evolution of the moments is analyzed. The control term, which is based on the mathematical thermostats, allows the derivation of nonlinear ordinary differential equations (Riccati-type equations) fulfilled by the moments.

The framework proposed here is certainly worth of future research activity concerning both its qualitative analysis and the application to modelling complex biological systems; see [7] and the references cited therein. Mathematical control in integrodifferential equations is already presented in the pertinent literatur; see, among others, [24].

The contents of the present paper are divided into three sections which follow this introduction. Section 2 highlights the essential hybrid mathematical settings. Section 3 deals with the analysis of the moments evolution. Finally, Section 4 is meant to future research perspectives.

2. The Underlying Hybrid Kinetic Framework

This section is meant to a concise description of an hybrid kinetic for active particle framework that constitutes a new paradigm for the modelling of complex multicellular systems characterized by nonconservative particle interactions. Specifically, these systems are composed by cells (active particles), with different genotypes and phenotypes, which are able to perform a task; for example, the immune system is constituted by cells that protect the human body from pathogens, and the connective tissue is made up of fibroblast cells that synthesize the extracellular matrix and collagen.

The mathematical structure consists of autonomous partial integrodifferential equations system with quadratic nonlinearity constituted by three different operators: the operator which models conservative interactions, the operator that takes care of nonconservative interactions, and the operator that ensures the control of the evolution of the moments. Each equation is a kinetic-type equation whose solution is the distribution function f(t, x, v, u) of cells that at time t ∈ [0, ∞[, possess the triplet of microscopic variables (x, v, u), where x ∈ D ⊂ ℝ3 is the space variable, v ∈ D ⊂ ℝ3 is the velocity variable, and u ∈ D ⊂ ℝ is the task-variable (activity). The elementary product f(t, x, v, u) d x  d v  du is the number of cells whose state, at time t, is in the elementary volume of the space of microscopic states dΩ = [x, x + d x]×[v, v + d v]×[u, u + du]. Let Ω = D × D × D be the domain of the all possible microscopic states and the Lebesgue measure on Ω. If f ∈ L 1[Ω, dμ ], then the following expectation 𝔼 0[f] of f: represents the global number, at time t, of the cells.

Higher order moments of f can be defined under suitable integrability assumption on f. Let ω = |x| | v| | u|, where |·| denotes the Euclidean distance of · and m, n, p ∈ ℕ. If, for instance, ω f ∈ L 1[Ω, dμ ], the (m, n, p)th-order moment of f is defined as follows: It is worth stressing that if f represents the joint distribution function, marginal distribution functions of f refer either to the distribution over the mechanical state or to distribution over the microscopic activity. These marginal distribution functions define the local quantities. For instance, the local quadratic activity moment at time t in x is computed as follows: and represent the energy expressed by the activity variable.

The task dynamics performed by cells is modified, at the time t, by binary interactions, which occur at the microscopic level and refer to the mutual actions between the cell with microscopic state (x, v, u) (test cell) and the cell with microscopic state (x*, v*, u*) (field cell), or the test cell and the cell with microscopic state (x ∗, v ∗, u ∗) (candidate cell). The candidate cell is the cell that reaches the microscopic state of the test cell after the interaction with the field cell.

2.1. The Hybrid Kinetic Setting at Equilibrium

Looking at the motion of cells, it is nonrestrictive to assume that velocity variables can attain discrete values v ∈ {v 1, v 2,…, v } and the activity variable is continuous u ∈ D . Therefore, the cell microscopic state includes both discrete and continuous variables. Moreover, we assume space homogeneity. Accordingly, the distribution function f of the system reads The mathematical framework thus consists of a set of evolution equations for f (t, u), which represents the distribution function of the test cell. Let f (t, u ∗) be the distribution function of the candidate cell and f (t, u*) the distribution function of the field cell.

The time evolution of the distribution function of test cell f (t, u) is modified when the candidate cell, with task-variable u ∗, after the interaction with the field cell, with task-variable u ∗, can acquire the task-variable u. The interactions in multicellular systems are conservative (modify the magnitude activity variable of the cells) and nonconservative (that modify the number of cells as a consequence of proliferation and destruction events occurring for natural birth/death and competition). Since f (t, u)du denotes the number of cells, at time t, such that u ∈ [u, u + du], then is a nonlinear product that refers to the interaction, at time t, between the number of candidate cells, with velocity v such that u ∗ ∈ [u ∗, u ∗ + du ∗] and the number of field cells, with velocity v , such that u* ∈ [u*, u* + du*]; the possibility of this interaction can be measured by introducing the nonnegative function ι (u ∗, u*) which represents the interaction rate; the probability that after this interaction the candidate cell undergoes a change in its microscopic state (that of test particle) is measured by the nonnegative function 𝒜 (u ∗, u*, u), which is a probability density with respect to u and then Bearing all above in mind, the (infinitesimal) result of the interaction reads and summing up with respect to all the candidate and field cells, we obtain the following operator which models the gain of test cell: Similarly, the loss of test cells is modeled by the following operator: Finally, if π (u, u*) is the net birth/death rate of the test cell due to the interaction with the field cell, then the operator which models nonconservative interactions reads The evolution equation of the distribution function f over the microscopic state can be derived by a balance equation of the inlet and outlet flows in the elementary volume [u, u + du] of the space of the microscopic states. The hybrid kinetic for active particles framework thus reads:

Definition 1

Let ι (u 1, u 2) : D × D → ℝ+, for i, j ∈ {1,2,…, n}, be the interaction rate between the u 1-cell distributed according to f (t, u 1) and the u 2-cell distributed according to f 2(t, u 2). Let 𝒜 (u 1, u 2, u) : D × D × D → ℝ+ be the probability density satisfying the property (7). A function f = f (t, u):(0, ∞) × D → ℝ+ is said to be the solution of (12) if

f (t, u) ∈ C((0, ∞), L 1(D ));

f is differentiable with respect to the variable t;

ι (u 1, u 2)𝒜 (u 1, u 2, u)f (t, u 1)f (t, u 2) is an integrable function with respect to the elementary measure du 1 du 2;

ι (u 1, u 2)f (t, u 2) is an integrable function with respect to the elementary measure du 2;

ι (u 1, u 2)π (u 1, u 2)f (t, u 2) is an integrable function with respect to the elementary measure du 2;

f satisfies (12) for all (t, u)∈(0, ∞) × D .

Setting f = (f 1(t, u), f 2(t, u),…, f (t, u)) ∈ ℝ, the (p, q)-order moment of the distribution function f(t, v, u), for p, q ∈ ℕ, is written as follows: In particular, the zero-order 𝔼 0,0 (density or mass), first-order 𝔼 1,1 (mean activation or linear momentum), and second-order 𝔼 2,2 (activation energy or kinetic energy) moments fulfill an important role depending on the system under consideration.

2.2. The Controlled Hybrid Kinetic Setting at Nonequilibrium

The mathematical framework (12) is concerned with multicellular systems at equilibrium. Nonequilibrium conditions occur when the system is subjected to external fields F (u) : D → ℝ+ at macroscopic scale. In this case, the kinetic framework reads The external field does work on the system thereby moving it away from equilibrium. Therefore, it follows uncontrolled increases of the activation energy (kinetic energy). The action of the Gaussian isokinetic thermostat is modeled, according to the Gauss' principle of least constrain [20, 25], by the following operator: which is a damping operator (thermostat operator) that is adjusted so as to control the activation energy. The introduction of the thermostat operator modifies the mathematical framework as follows: where 𝒞 [f , f ](t, u) = 𝒢 [f , f ](t, u) − ℒ [f , f ](t, u) is the operator for the conservative interactions. In what follows, we refer to framework (16) as the controlled kinetic framework with conservative and nonconservative interactions.

Definition 2

Let F = F (u), u ∈ D , be an external force field differentiable with respect to u; ι (u 1, u 2) : D × D → ℝ+, for i, j ∈ {1,2,…, n}, interaction rate between the u 1-cell distributed according to f (t, u 1) and the u 2-cell distributed according to f 2(t, u 2); consider 𝒜 (u 1, u 2, u) : D × D × D → ℝ+ to be the probability density satisfying the property (7). A function f = f (t, u):(0, ∞) × D → ℝ+ is said to be the solution of the model (16) if

f (t, u) ∈ C((0, ∞), L 1(D ));

f is differentiable with respect to the variables t and u;

uf is an integrable function with respect to the elementary measure du;

ι (u 1, u 2)𝒜 (u 1, u 2, u)f (t, u 1)f (t, u 2) is an integrable function with respect to the elementary measure du 1 du 2;

ι (u 1, u 2)f (t, u 2) is an integrable function with respect to the elementary measure du 2;

ι (u 1, u 2)π (u 1, u 2)f (t, u 2) is an integrable function with respect to the elementary measure du 2;

𝒯 [F , f] is differentiable with respect to the variable u;

f satisfies (16) for all (t, u)∈(0, ∞) × D .

Remark 3

The theorem of existence and uniqueness of the solution for the controlled kinetic framework (16) has been obtained in [7] when the nonconservative operator 𝒩 is equal to zero (conservative interactions only). The proof of the theorem can be adapted in order to obtain existence and uniqueness of the solution also for the nonconservative interactions case. Nevertheless, global existence may not occur. This is a work in progress and results will be reported in due course.

The depicted hybrid controlled kinetic framework (16) is quite general and can be exploited to originate specific models for multicellular systems by acting on the specific forms of the grid velocity, interaction rate ι , the probability density 𝒜 , the net rate of birth/death π , and the external force F .

3. Differential Equations for the Moments

This section is concerned with the derivation of differential equations for the moments. Let 𝔼 1[f] be the following moment: Let μ(t) be the following function: where

The following result holds true.

Theorem 4

Let ι = ι, π = π, and F = F be real constants. If there exists a nonnegative solution f of the controlled kinetic framework (16) such that f (t, u) = 0 as u ∈ ∂D , then the 1th-order moment 𝔼 1[f](t) is solution of the following Riccati nonlinear ordinary differential equation:

Proof

The interaction operator 𝒥[f , f ] can be written as follows: Multiplying both sides of 𝒥 [f , f ] by u and integrating over D , we have Summing with respect to j, multiplying by v , and summing with respect to i, we obtain Multiplying by u and v the second term of the left hand side of (16), integrating with respect to the activity variable, performing integration by parts and summing with respect to i, we have and then the proof.

According to Theorem 4, the solution of the Riccati equation (21) can be obtained as follows. The Riccati equation reads Setting if is a solution of (26), the general integral can be written as where λ(t) is solution of A nonnegative and constant solution of (26) is Therefore, the solution of (26) can be written as follows: where

The next theorem gives the evolution equation for all moments where p is an odd number.

Theorem 5

Let p ∈ ℕ, q be an odd number and t ≥ 0. Then, the (p, q)th-order moment of the distribution function f satisfies the following Riccati nonlinear ordinary differential equation: where Moreover, if   𝔼 [f](t) is initially bounded, it remains bounded for all t > 0.

The proof follows by multiplying both sides of (14) by u and performing integration by parts on the control term.

4. Research Perspectives

The controlled kinetic framework proposed in this paper allows the derivation of specific models for multicellular systems characterized by nonconservative interactions. This framework belongs to the class of thermostated kinetic for active particles models.

The mathematical framework (16) can be further generalized in order to include the role of mutations; see Nowak [26]. This is an important issue in the cancer modeling [27, 28].

A future research perspective is the generalization of the mathematical framework (16) to open systems subjected to external actions at the microscopic scale, for example, the role that the outer environment has in the whole dynamics [29].

Perspectives include also the introduction of stochastic terms that model jump processes in the activity or in the velocity variable; see paper [30] and the references section.

Moreover, the proof of the existence of solutions to the stationary problem is missing. This proof has been gained for the conservative interactions case; see [31].

An important research perspective is the development of suitable asymptotic limits for deriving macroscopic equations for the evolution of moments; see [17] and the references section. The derivation of these equations is based on suitable assumptions on the operators in order to obtain convergence results. Optimization investigation can also be performed on the solutions of these macroscopic equations, such as regularization and embedding results [32].

It is worth stressing that the analysis of the moments performed in this paper can be straightforwardly applied also for the thermostated framework proposed in [9].