Hybrid models in biomedical applications

КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ 2019 Т. 11 № 2 С. 287–309 DOI: 10.20537/2076-7633-2019-11-2-287-309 АНАЛИЗ И МОДЕЛИРОВАНИЕ СЛОЖНЫХ ЖИВЫХ СИСТЕМ УДК: 519.8 Гибридные модели в биомедицинских приложениях Н. М. Бессонов1,a , Г. А. Бочаров2,b , А. Бушнита3,c , В. А. Вольперт4,5,6,7,d 1 Институт проблем машиноведения РАН, Россия, 199178, г. Санкт Петербург, Большой проспект В.О., д. 61 2 Институт вычислительной математики им. Г. И. Марчука РАН, Россия, 119333, г. Москва, ул. Губкина, д. 8 3 Факультет информационных технологий, Уппсальский Университет, Швеция, SE-752 37, г. Уппсала, ул. Лагерхидсваген, д. 2 4 Институт Камиля Жордана, УМР 5208 НЦНИ, Университит Лион 1, Франция, г. Вилеурбан, бульвар 11 ноября 1918 года, д. 43 5 ИНРИА, Институт Камиля Жордана, Университит Лион 1, Франция, 69603, г. Вилеурбан, бульвар Нильса Бора, д. 56 6 Российский университет дружбы народов, Россия, 117198, г. Москва, ул. Миклухо-Маклая, д. 6 7 Центр Понселе, УМИ 2615 НЦНИ, Россия, 119002, г. Москва, Большой Влассиевский пер., д. 11 E-mail: a nickbessonov@yahoo.com , b bocharov@m.inm.ras.ru, c anass.bouchnita@it.uu.se, d volpert@math.univ-lyon1.fr Получено 28.09.2018, после доработки — 06.11.2018. Принято к публикации 24.01.2019. В статье представлен обзор недавних работ по гибридным дискретно-непрерывным моделям в динамике клеточных популяций. В этих моделях, широко используемых в биологическом моделировании, клетки рассматриваются как отдельные объекты, которые могут делиться, умирать, дифференцироваться и двигаться под воздействием внешних сил. В простейшем представлении клетки рассматриваются как мягкие сферы, их движение описывается вторым законом Ньютона для их центров. В более полном представлении могут учитываться геометрия и структура клеток. Судьба клеток определяется концентрациями внутриклеточных веществ и различных веществ во внеклеточном матриксе, таких как питательные вещества, гормоны, факторы роста. Внутриклеточные регуляторные сети описываются обыкновенными дифференциальными уравнениями, а внеклеточные концентрации — уравнениями в частных производных. Мы проиллюстрируем применение этого подхода некоторыми примерами, в том числе бактериальными филаметами и ростом раковой опухоли. Далее будут приведены более детальные исследования эритропоэза и иммунного ответа. Эритроциты произодятся в костном мозге в небольших структурах, называемых эритробластными островками. Каждый островок образован центральным макрофагом, окруженным эритроидными предшественниками на разных стадиях зрелости. Их выбор между самообновлением, дифференцировкой и апоптозом определяется регуляцией ERK/Fas и фактором роста, производимым макрофагами. Нормальное функционирование эритропоэза может быть нарушено развитием множественной миеломы, злокачественного заболевания крови, которое приводит к разрушению эритробластических островков и к развитию анемии. Последняя часть работы посвящена применению гибридных моделей для изучения иммунного ответа и развития вирусной инфекции. Представлена двухмасштабная модель, включающая лимфатический узел и другие ткани организма, включая кровеносную систему. Ключевые слова: клеточные полуляции, дискретно-непрерывные модели, эритропoэз, иммунный ответ Работа частично поддержана программой «Университет РУДН 5-100». Исследования по гибридным моделям иммунного ответа (секция 4) поддержаны Российским научным фондом (грант № 18-11-00171). c 2019 Николай Михайлович Бессонов, Геннадий Алексеевич Бочаров, Анаcс Бушнита, Виталий Айзикович Вольперт  Статья доступна по лицензии Creative Commons Attribution-NoDerivs 3.0 Unported License. Чтобы получить текст лицензии, посетите вебсайт http://creativecommons.org/licenses/by-nd/3.0/ или отправьте письмо в Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. COMPUTER RESEARCH AND MODELING 2019 VOL. 11 NO. 2 P. 287–309 DOI: 10.20537/2076-7633-2019-11-2-287-309 ANALYSIS AND MODELING OF COMPLEX LIVING SYSTEMS UDC: 519.8 Hybrid models in biomedical applications N. M. Bessonov1,a , G. A. Bocharov2,b , A. Bouchnita3,c , V. A. Volpert4,5,6,7,d 1 Institute of Problems of Mechanical Engineering, Russian Academy of Sciences, V.O., 61 Bolshoj pr., St. Petersburg, 199178, Russia 2 Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences, 8 Gubkina st., Moscow, 119333, Russia 3 Department of Information Technology, Uppsala University Lägerhyddsvägen 2, Uppsala, SE-752 37, Sweden 4 Institut Camille Jordan, UMR 5208 CNRS, Université Lyon 1, 43 Boulevard du 11 Novembre 1918, Villeurbanne, 69622, France 5 INRIA, Institut Camille Jordan, Université Lyon 1, 56 Boulevard Niels Bohr, Villeurbanne, 69603, France 6 Peoples’ Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya st., Moscow, 117198, Russia 7 Poncelet Center, UMI 2615 CNRS, 11 Bolshoi Vlassievskii, Moscow, 119002, Russia E-mail: a nickbessonov@yahoo.com , b bocharov@m.inm.ras.ru, c anass.bouchnita@it.uu.se, d volpert@math.univ-lyon1.fr Received 28.09.2018, after completion — 06.11.2018. Accepted for publication 24.01.2019. The paper presents a review of recent developments of hybrid discrete-continuous models in cell population dynamics. Such models are widely used in the biological modelling. Cells are considered as individual objects which can divide, die by apoptosis, differentiate and move under external forces. In the simplest representation cells are considered as soft spheres, and their motion is described by Newton’s second law for their centers. In a more complete representation, cell geometry and structure can be taken into account. Cell fate is determined by concentrations of intra-cellular substances and by various substances in the extracellular matrix, such as nutrients, hormones, growth factors. Intra-cellular regulatory networks are described by ordinary differential equations while extracellular species by partial differential equations. We illustrate the application of this approach with some examples including bacteria filament and tumor growth. These examples are followed by more detailed studies of erythropoiesis and immune response. Erythrocytes are produced in the bone marrow in small cellular units called erythroblastic islands. Each island is formed by a central macrophage surrounded by erythroid progenitors in different stages of maturity. Their choice between self-renewal, differentiation and apoptosis is determined by the ERK/Fas regulation and by a growth factor produced by the macrophage. Normal functioning of erythropoiesis can be compromised by the development of multiple myeloma, a malignant blood disorder which leads to a destruction of erythroblastic islands and to sever anemia. The last part of the work is devoted to the applications of hybrid models to study immune response and the development of viral infection. A two-scale model describing processes in a lymph node and other organs including the blood compartment is presented. Keywords: cell populations, discrete-continuous models, erythropoiesis, immune response Citation: Computer Research and Modeling, 2019, vol. 11, no. 2, pp. 287–309 (Russian). The work was partially supported by the “RUDN University Program 5-100”. The research on hybrid modelling of immune response (Section 4) was supported by the Russian Science Foundation (grant no. 18-11-00171). c 2019 Nikolai M. Bessonov, Gennady A. Bocharov, Anass Bouchnita, Vitaly A. Volpert  This work is licensed under the Creative Commons Attribution-NoDerivs 3.0 Unported License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nd/3.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. Hybrid models in biomedical applications 289 1. Discrete-continuous models of cell populations This review is devoted to the modelling of cell population dynamics with so-called hybrid models where cells are considered as discrete or individual based objects while intra-cellular and extra-cellular biochemical substances are described with continuous models. It is in this sense that we understand here hybrid models. Hybrid models can be based on cellular automata or various lattice or off-lattice approaches (see [Anderson, 2007; Anderson et al., 2007a; Anderson et al., 2007b; Drasdo, 2007] and the references therein). Cells can interact with each other and with the surrounding medium mechanically and biochemically, they can divide, differentiate and die due to apoptosis. Cell behavior is determined by intra-cellular regulatory networks described by ordinary differential equations and by extra-cellular bio-chemical substances described by partial differential equations. In order to describe mechanical interaction between cells, we restrict ourselves here to the simplest model where cells are represented as elastic balls. Consider two elastic balls with the centers at the points x1 and x2 and with the radii, respectively, r1 and r2 . If the distance d12 between the centers is less then the sum of the radii, r1 + r2 , then there is a repulsive force between them f12 which depends on the distance d12 . If a particle with the center at xi is surrounded by several other particles with the centers at the points x j , j = 1, . . . , k, then we consider the pairwise forces fi j assuming that they are independent of each other. This assumption corresponds to small deformation of the particles. Hence,  we find the total force Fi acting on the i-th particle from all other particles, Fi = ji fi j . The motion of the particles can now be described as the motion of their centers by Newton’s second law  f (di j ) = 0, (1.1) m ẍi + µm ẋi − ji where m is the mass of the particle, the second term in the left-hand side describes the friction by the surrounding medium. Dissipative forces can also be written in a different form. This is related to dissipative particle dynamics [Karttunen et al., 2004]. More complex models of motion can be considered if the geometry of the moving object (particle, cell, organism) and its properties should be taken into account (see, e.g., [Bessonov et al., 2014; Tosenberger et al., 2016; Bessonov et al., 2015]). Intra-cellular regulatory networks for the i-th cell are described by a system of ordinary differential equations dui = F(ui , u), (1.2) dt where ui is a vector of intra-cellular concentrations, u is a vector of extra-cellular concentrations, F is the vector of reaction rates which should be specified for each particular application. The concentrations of the species in the extra-cellular matrix is described by the diffusion equation ∂u = D ∆u + G(u, c), ∂t (1.3) where c is the local cell density, G is the rate of consumption or production of these substances by cells. These species can be either nutrients coming from outside and consumed by cells or some other bio-chemical products consumed or produced by cells. In particular, these can be hormones or other signaling molecules that can influence intra-cellular regulatory networks. In some cases, convective motion of the medium should be taken into account. There are various details of this model related to cell division, the force fi j and cell displacement, relation between discrete and continuous models of cell population, more complex cell geometry, influence of stochasticity and so on (see [Bessonov et al., 2006]– [Bessonov et al., 2010a] for more details). In the next section we present some model examples. Section 3 is devoted to erythropoiesis and multiple myeloma, in Section 4 we review hybrid models of immune response. 2019, Т. 11, № 2, С. 287–309 290 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert 2. Model examples 2.1. Self-renewal and apoptosis We begin with the 1D model example where cells can move along the line. The coordinates xi in equations (1.1) are real numbers. Each cell can divide or die by apoptosis. After division a cell gives two cells identical to itself. We suppose that cell division and death are influenced by some bio-chemical substances produced by the cells themselves. We consider the case where there are two such substances, u and v. We come to the following system of equations: ⎧ ⎪ ∂2 u du ⎪ ⎪ ⎪ + b1 c − q1 u, = d ⎪ 1 ⎪ ⎨ dt ∂x2 (2.1) ⎪ ⎪ ⎪ ∂2 v dv ⎪ ⎪ ⎪ = d2 2 + b2 c − q2 v. ⎩ dt ∂x These equations describe the evolution of the extracellular concentrations u and v with their diffusion, production terms proportional to the concentration of cells c and with the degradation terms. We note that cells are considered here as point sources with a given rate of production of u and v. The cell concentration is understood as a number of such sources in a unit volume. In numerical simulations, where cells have a finite size, we consider them as distributed sources and specify the production rate for each node of the numerical mesh. Intra-cellular concentrations ui and vi in the i-th cell are described by the equations: ⎧ du ⎪ i ⎪ ⎪ = k1(1) u(x, t) − k2(1) ui (t) + H1 , ⎪ ⎪ ⎨ dt (2.2) ⎪ ⎪ ⎪ dvi ⎪ ⎪ ⎩ = k1(2) v(x, t) − k2(2) vi (t) + H2 . dt Here and in what follows we write equations for intra-cellular concentrations neglecting the change of the cell volume. This approximation is justified since the volume changes only twice before cell division and this change is relatively slow. The first term in the right-hand side of the first equation shows that the intra-cellular concentration ui grows proportionally to the value of the extra-cellular concentration u(x, t) at the space point xi where the cell is located. It is similar for the second equation. These equations contain the degradation terms and constant production terms, H1 and H2 . When a new cell appears, we put the concentrations ui and vi equal zero. If the concentration ui attains some critical value uc , then the cell divides. If vi becomes equal vc , the cell dies. Consider first the case where k1(1) = k1(2) = k2(1) = k2(1) = 0. Then ui and vi are linear functions of time which reach their critical values at some times t = τu and t = τv , respectively. If τu < τv , then all cells will divide with a given frequency, if the inequality is opposite, then all cells will die. Next, consider the case where k1(1) is different from zero. If it is positive, then cells stimulate proliferation of the surrounding cells, if it is negative, they suppress it. Both cases can be observed experimentally. We restrict ourselves here by the example of negative k1(1) . All other coefficients remain zero. Therefore, cells have a fixed life time τv . If they do not divide during this time, they die. We carry out the 1D simulation where cells can move along the straight line. Initially, there are two cells in the middle of the interval. Figure 1 shows the evolution of this population in time. For each moment of time (vertical axis) we have the positions of cells (horizontal axis) indicated with blue points. The evolution of the cell population in Figure 1 (top, left) can be characterized by two main properties. First of all, it expands to the left and to the right with approximately constant speed. Second, the total population consists of relatively small sub-populations. Each of them starts from КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 291 Figure 1. Dynamics of cell population in the case where cells either self-renew or die by apoptosis. Cells are shown with blue dots. Horizontal axis shows cell position, vertical axis — time a small number of cells. Usually, these are two cells at the right and at the left of the previous subpopulation. During some time, the sub-population grows, reaches certain size and disappears giving birth to new sub-populations. This behavior can be explained as follows. The characteristic time of cell division is less than of cell death. When the cell sub-population is small, the quantity of u is also small, and its influence on cell division is not significant. When the sub-population becomes larger, it slows down cell division because of growth of u. As a result the sub-population disappears. The outer cells can survive because the level of u there is less. The geometrical pattern of cell distribution for these values of parameters reminds Serpinsky carpet (Figure 1 top, left), an example of fractal sets. The pattern of cell distribution depends on the parameters. Other examples are shown in Figure 1. The cell populations in Figure 1 (bottom) remain bounded, and the patterns are time periodic. The simulations presented here do not use the extra-cellular variable v. Instead of the variable u, which decelerate cell proliferation, we can consider v assuming that it accelerates cell apoptosis. In this case, qualitative behavior of cell population is similar. 2.2. Bacteria filaments The model examples considered in the previous section should be adapted in the case of specific biological applications. In this section we will consider one of the simplest examples of morphogenesis, filaments of bacteria with self-renewing and differentiated cells. Filament of bacteria represents a chain of cells connected to each other. They can divide producing identical cells. When the cells in the filament lack nitrogen, some of them differentiate, some other remain in their original form. Differentiated cells do not divide. They are located periodically in the filament being separated by a given number of undifferentiated cells. In the case of anabaena filament, the intra-cellular regulation which determines cell differentiation is shown in Figure 2 [Sakr et al., 2006]. One of the earliest steps of heterocyst differentiation is the accumulation of 2-oxoglutarate (2-OG), which constitutes a signal of nitrogen starvation. It initiates production of the protein HetR which plays a key role in this regulation. First of all, its amplifies its own production. Next, it initiates production of the protein PatS. Finally, when the concentration of HetR becomes sufficiently high, the cell differentiates into heterocyst. On the other hand, PatS suppresses HetR and, as a consequence, cell differentiation. Moreover, PatS can diffuse between the neighboring cells. This competition between HetR and PatS is supposed to determine the pattern of differentiated cells. We emphasize that if a cell differentiates, then PatS prevents differentiation of the neighboring cells but not of the cell itself. So the proposed mechanism should capture this property of the system. 2019, Т. 11, № 2, С. 287–309 292 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert The intra-cellular concentrations are described by the following system of equations: ⎧ dui ⎪ ⎪ ⎪ = Hu , ⎪ ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dGi ⎪ ⎪ ⎪ = Hg , ⎪ ⎨ dt ⎪ ⎪ ⎪ ⎪ dhi ⎪ ⎪ ⎪ = k1(h)Gi (t) + k2(h) h2i (t) − k3(h) pi (t)hi (t), ⎪ ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ dpi ⎪ ⎪ ⎩ = d(pi−1 − pi ) + d(pi+1 − pi ) + k1(p) hi (t) − k2(p) pi (t), dt (2.3) where ui denotes the intra-cellular concentration of the cell division protein, FtsZ, in the i-th cell, Gi is the concentration of 2-oxoglutarate, hi is the self-enhancing differentiation regulatory protein HetR, pi is the inhibitor encoded by the gene PatS, which is dependent on HetR for production. PatS is synthesized in the developing heterocyst and diffuses to the neighboring vegetative cells. Since the right-hand side in the equation for ui is constant, then ui = Hu t. If this concentration reaches a critical value uc , the cell divides. The time of cell proliferation is τ p = uc /Hu . During this time, cell grows linearly and becomes twice bigger when its age equals τ p . The daughter cells are twice as small as the mother cell. Figure 2. Anabaena growth. Differentiated cells are red, undifferentiated yellow. The black circle inside cells shows their incompressible part. Color version of the picture is available on the journal website Cell differentiation occurs if the concentration of hi equals some critical value hc . Therefore, the choice between cell division and differentiation depends on what critical value is reached first. The rate of production of hi depends on Gi , on hi itself with quadratic self-acceleration and on pi which downregulates it. The right-hand side of the equation for pi contains the flux terms from the surrounding cells which are similar to the discretized diffusion equation, the production term proportional to hi and the degradation term. After differentiation, the cell does not grow anymore, it cannot divide and it keeps a constant level of pi . We need to specify the initial values of the concentrations in a new cell. When a cell divides, the initial conditions in the daughter cells are as follows: ui (t0 ) = 0, Gi (t0 ) = G0 , hi (t0 ) = 0, pi (t0 ) = pm i , 2 (2.4) where pm i quantity of PatS in the mother cell. Hence the initial concentrations of FtsZ and HetR are zero, 2-OG takes on some constant given value, and PatS equals half of the concentration in the mother cell. The latter is important for the pattern of heterocyst differentiation. Indeed, if a cell is located sufficiently close to a differentiated cell, then the concentration of PatS in it is sufficiently high. When it divides, the concentration of PatS in the daughter cell, though two times less, can be sufficient to keep the level of HetR less than the critical level. If the mother cell is located far from differentiated cells, then the initial concentration of PatS in the daughter cell is close to zero, and it cannot prevent growth of HetR up to the critical value which determines differentiation. КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 293 An example of numerical simulations is shown in Figure 2. We begin with 5 cells, one differentiated cell in the center and two undifferentiated cells from each side. It is also possible to begin with all undifferentiated cells. In this case we need to introduce random perturbations in the production terms for HetR in order to give an advantage to one of the cells to differentiate. Otherwise, if all cells are initially identical, then they will remain identical. Hence, we need to introduce some asymmetry either in initial conditions or in the equations. Undifferentiated cells begin to grow, then they divide. For t = 92.8, the most left and the most right cells differentiate. Undifferentiated cells in between divide once again and after some time a cell at the center of the subpopulation of undifferentiated cells becomes differentiated. The process of growth, division and differentiation continues in the same manner. The number of undifferentiated cells between two nearest differentiated cells is not exactly constant because of the parity due to division. Let us explain this in more detail. At t = 92.8 there are 7 undifferentiated cells from each side. After division, there are 14 of them. A new differentiated cell splits them into two groups of 7 and 6 cells. After the next division, there are 14 and 12 cells. The next differentiation will produce the groups of 7, 6, 6, and 5 cells. Hence, after each cycle, there is a group of undifferentiated cells smaller than previously. At some point, there will be a group small enough such that it will not give a differentiated cell. Then the cells in this group will divide producing a group of cells twice larger. 2.3. Tumor growth In this section we consider the complete model (1.1)–(1.3) assuming that u is a scalar variable. It describes the concentration of nutrients which diffuse from the boundary of the domain and which are consumed by cells inside the domain. We consider the equation ∂u = D ∆u − kcu, ∂t (2.5) where c is cell concentration and k is a positive parameter. The rate of nutrient consumption is proportional to the product of the concentrations. The form of the domain and of the boundary conditions will be specified below. Next, we consider the scalar intra-cellular variable ui where the subscript i corresponds to the cell number. It is described by the equation dui = k1 u(xi , t) − k2 ui . dt (2.6) The first term in the right-hand side of this equation shows that the intra-cellular concentration ui grows proportionally to the value of the extra-cellular concentration u(x, t) at the space point xi where the cell is located. The second term describes consumption or destruction of ui inside the cell. We suppose that the cell radius ri grows proportionally to the increase of ui :  dui dri = max ,0 . (2.7) dt dt The initial value of the radius for each new cell is r0 , the maximal radius rm . When it is reached, the cell divides. If the cell does not divide before its maximal age, then it dies. The maximal cell age is a parameter of the problem. Consider a circular domain Ω and complete equation (2.5) by the boundary condition u = 1 at the boundary ∂Ω. We put a single cell in the center of the domain and begin numerical solution of system (2.5)–(2.7). The results of the simulations are shown in Figure 3 for several consecutive moments of time. The grey 2D surface shows the spatial distribution of the concentration u(x, t). In the beginning it equals 1 everywhere in the domain. The constants k1 and k2 , k1 > k2 are chosen in such 2019, Т. 11, № 2, С. 287–309 294 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert Figure 3. Consecutive moments of tumor growth. It starts with a single cell at the center of the circle. Living cells are shown in red, dead cells in black. 2D grey surface shows the level of nutrients. Color version of the picture is available on the journal website a way that the intra-cellular concentration ui growth. Consequently, the radius of the cell also grows and after some time the cell divides. The new cells also consume nutrients, grow and divide. The part of the domain filled by cells forms a disk while the concentration u(x, t) decreases in the center of the domain (second figure in the upper row). Hence, the right-hand side of equation (2.7) also decreases, the intra-cellular concentration stops growing or even decreases, and cells cannot divide before their maximal age τ. As a result, they die and form the black region in the center. Living cells shown in red form a narrow external layer. The region filled by cells grows in time and finally approaches to the boundary of the domain. Dynamics of the cell population can be more complex if we decrease the maximal life time τ. Cells now have less time to accumulate enough nutrients for division. In this case, even a small decrease in nutrient concentration can become crucial from the point of view of the choice between proliferation and apoptosis. In the beginning cells form, as before, a circular region with living cells outside and the core formed by dead cells. The layer of living cells is narrower than in the previous example. After that the domain loses its radial symmetry resulting in the appearance of spatial patterns of the growing tumor (not shown). 3. Erythropoiesis and multiple myeloma 3.1. Organization of erythropoiesis Erythropoiesis (production of red blood cells) represents continuous throughout life process, maintaining an optimal number of circulating red cells and tissue oxygen tension. Erythropoiesis happens in the bone marrow where erythroid progenitors, immature blood cells with abilities to proliferate and differentiate, undergo a series of maturation stages to produce erythroblasts (mature progenitors) and then reticulocytes which subsequently enter the bloodstream and mature into red КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 295 cells (erythrocytes). At every step of this differentiation process, erythroid cells can die by apoptosis (programmed cell death) or self-renew [Bauer et al., 1999; Gandrillon, 2002; Gandrillon et al., 1999; Pain et al., 1991]. Numerous external regulations control cell fate by modifying the activity of intracellular proteins. Erythropoietin (Epo) is the hormone synthesized in the kidney in response of decreased in tissue oxygen tension. Epo promotes survival of early erythroblast subsets by negative regulation of their apoptosis through the action on the death receptor Fas [Koury, Bondurant, 1990]. Moreover, Epo induces self-renewal [Rubiolo et al., 2006; Spivak et al., 1991] in cooperation with glucocorticoids [Bauer et al., 1999; Gandrillon, 2002; Gandrillon et al., 1999; Pain et al., 1991] and some intracellular autocrine loops [Gandrillon et al., 1999; Sawyer, Jacobs-Helber, 2000]. It has been shown [Rubiolo et al., 2006], that intracellular proteins Erk and Fas are involved in two antagonist loops. Previously considered by the authors [Crauste et al., 2010], Erk (from the MAPK family) inhibits Fas (a TNF family member) and self-activates inducing cell proliferation depending on Epo levels, whereas Fas inhibits Erk inducing apoptosis and differentiation. Fate of cell depends on the level of each protein. In addition to global feedbacks, there is local feedback control through cell-cell interaction, during which reticulocyte-produced Fas-ligand binds to the membrane protein Fas inducing both differentiation and death by apoptosis [De Maria et al., 1999]. All the process of erythroid maturation occurs in the context of erythroblastic islands, the specialized niches of bone marrow, in which erythroblasts surround a central macrophage which appears to influence their proliferation and differentiation [Chasis, Mohandas, 2008; Tsiftsoglou et al., 2009]. For years, however, erythropoiesis has been studied under the influence of erythropoietin, which may induce differentiation and proliferation in vitro without the presence of the macrophage. Hence, the roles of the macrophage and the erythroblastic island have been more or less neglected. In the following paper proposed model was validated by comparing with experimental data during stress erythropoiesis as hypoxia and different forms of induced anaemia [Wichmann et al., 1989]. In 1995, Bélair et al. proposed age-structured model of erythropoiesis where erythropoietin causes differentiation, without taking into account erythropoietin control of apoptosis found out in 1990 by Koury and Bondurant in [Koury, Bondurant, 1990]. In 1998 Mahaffy et al. [Mahaffy et al., 1998] expanded this model by including the apoptosis possibility. Age-stuctured model is detailed in [Ackleh et al., 2006] with assumption that decay rate of erythropoietin depends on the number of precursor cells. In [Crauste et al., 2008] Crauste et al. included in the model the influence of Epo upon progenitors apoptosis and showed the importance of erythroid progenitor self-renewing by confronting their model with experimental data on anaemia in mice. Multiscale approache, that include both erythroid progenitor dynamics [Crauste et al., 2010; Crauste et al., 2008] and intracellular regulatory networks dynamics, give insight in the mechanisms involved in erythropoiesis. We present here a hybrid discrete-continuous modelling approach [Bessonov et al., 2011; Bessonov, 2010b] to bring together intracellular and extracellular levels as well as spatial aspects, in order to study the importance of spatial structure of erythroblastic island in erythropoiesis. In previous models of erythropoiesis we considered the process in human bone marrow erythroblastic islands [Fischer et al., 2012; Kurbatova et al., 2011] with local feedback, mediated by reticulocytes exhibiting a Fas-based cytotoxicity against the immature erythroblasts via Fas-ligand induction [De Maria et al., 1999]. A more complete model is presented here with the production of red blood cell in mice and taking into account the aspect of co-expressing of Fas and Fas-ligand by immature erythroid progenitors, particularly in the spleen [Koulnis et al., 2011]. 3.2. Model of erythroblastic islands Some erythroid cells produce Fas-Ligand, F L , which influences the surrounding cells by increasing intracellular Fas activity. These are immature cells in murine erythropoiesis [Koulnis et al., 2011] and more mature cells in human erythropoiesis [De Maria et al., 1999]. On the other 2019, Т. 11, № 2, С. 287–309 296 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert hand, macrophages produce a growth factor G, which stimulates proliferation. The concentrations of Fas-ligand and of the growth factor in the extracellular matrix are described by the reaction-diffusion equations: ∂F L = D1 ∆F + W1 − σ1 F L , ∂t ∂G = D2 ∆G + W2 − σ2G. ∂t Here W1 and W2 are the rates of production of the corresponding factors, the last terms in the right-hand sides of these equations describe their degradation. Different systems of ordinary differential equations can be used for the intracellular regulation. The action of Fas-ligand and of the growth factor on the intracellular regulation is modeled by the value parameters. They become linear functions of quantities of this two substances: γ1 = γ10 + γ11 B, γ2 = γ20 + γ21 E, γ3 = γ30 + γ31 F L with B — BMP4 produced by macrophage, E — Epo, F L — Fas-ligand. Simplified model of intracellular regulation in erythroid progenitors is described by ordinary differential equations. Behavior of cells depends on the values of u, v, w. At the end of a cellular time cycle, if u < v, progenitors become a differentiated cell, else if u > v, progenitors choose self-renewal. If w > wcrit at any point in the life of a progenitor, the cell die by apoptosis. du = γ1 , dt dv = γ2 (1 − β1 uv), dt dw = γ3 (1 − β2 uw). dt (3.1) Figure 4. Two stages of the development of multiple myeloma. Myeloma cells form multiple tumors and destroy erythroblastic islands mechanically and biochemically. In the beginning of myeloma development, erythroblastic islands are weakly influenced and continue their functioning (upper figure). After some time, tumor fills an important part of the marrow and all islands are destroyed (lower figure). An erythroblastic island is formed by a large cell, the macrophage produces GF (green halo), immature cells (progenitors) in yellow, mature cells (differentiated cells) in blue. Reprinted from [Bouchnita et al., 2016]. Color version of the picture is available on the journal website Example of numerical simulations is shown in Figure 4. Macrophage, large cell in the center of the island produces a growth factor that influences proliferation of erythroid progenitors around it. These cells can differentiate into reticulocytes or dies by apoptosis. Reticulocytes produce Fas-ligand upregulating apoptosis of erythroid progenitors. КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 297 3.3. Multiple myeloma The development of (multiple myeloma) MM leads to several harmful clinical conditions including anemia, renal failure, recurrent infections, hypercalcemia, and osteoporosis with bone fractures. These pathological conditions can frequently result in the death of the patient. Multiple myeloma tumor destroys erythroblastic islands undermining normal functioning of erythropoiesis (Fig. 4, right). The evolution of myeloma cells from the clonal premalignant plasma cells due to monoclonal gammopathy of unknown significance (MGUS) through the intermediate stage of smoldering myeloma to the malignant MM stage is mediated by multiple sequential genetic changes, including chromosomal translocations, hyperdiploidy, and mutations, which permit independent growth and spread of MM in the bone marrow [Palumbo, Anderson, 2011; Morgan et al., 2012]. A wide variety of genetic changes involving many different genes have been documented in MM cases. In addition to these sequential genetic changes, MM cells require specific interactions with the non-hematopoietic cells of the bone marrow including the stromal cells (BMSCs), osteoblasts, osteoclasts, and cells associated with vascular supply of the marrow [Palumbo, Anderson, 2011; Morgan et al., 2012]. These marrow microenvironment interactions with the MM cells include direct binding of cell surface adhesion proteins and their binding partners on the other cell types as well as diffusible, soluble molecules that are secreted by one cell type and are bound and internalized by another cell type. These soluble molecules are often chemokines and cytokines produced by other marrow cells, and their receptors are expressed on the MM cells. However, specific molecules produced by MM cells indirectly affect growth of the malignancy by inducing localized resorption of bone, death of hematopoietic cells, and expansion of vessels supplying blood to the MM. The combined effects of sequential genetic changes within the evolving myeloma cells and interactions of the MM cells with the marrow microenvironment affect the migration, proliferation, and survival of the MM cells. The model is based on the direct effects of sequential genetic changes and marrow microenvironmental chemokine and cytokine activity that influence the chemotaxis, proliferation and survival of MM, but does not include the MM effects on bone resorption, hematopoietic cell loss, or development of specialized vasculature. The complex multiple genetic changes in MM cells and the numerous cell–cell and cytokine-mediated interactions between myeloma cells and their marrow microenviroment are simplified in the model so that four related but evolving clones develop in a process termed intra-clonal heterogeneity [Brioli et al., 2014; Melchor et al., 2014] (Figure 5). Competition among these four MM clones is based on differences in cellular growth and survival rates and interactions with the marrow microenvironment. This competition results in predominance of the more fit clones and decline and ultimate extinction of the less fit ones. 4. Immune response The immune system provides the defense of a host organism against foreign pathogens and tumor development, and plays an active role in tissue and organ regeneration. Deviations of the system functioning from normal physiological activity lead to diseases with various pathologies including autoimmune- and cancer processes. The modern era of research in immunology is characterized by an unprecedented level of detail about its numerous components functioning together as a whole network-type system. However, pure empirical analyses of the behavior of the immune system, and of its response to external perturbations, are limited to a static description of its components and the connections between them. There is a demand for the development of high-resolution detailed mathematical models and their integration into experimental and clinical research to provide a mechanistic tool for the description, analysis and prediction of immune process dynamics under specified conditions. (The material of the section is based on our previous studies [Bouchnita et al., 2017b; Bouchnita et al., 2017c].) 2019, Т. 11, № 2, С. 287–309 298 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert Figure 5. (a) The intracellular regulation of myeloma cells as described in the model. Bone marrow stromal cells (BMSCs) secrete the cytokines insulin-like growth factor 1 (IGF-1), interleukin-6 (IL-6), and chemokine stromal cell-derived factor 1 (SDF-1) which are necessary to the survival, homing and proliferation of myeloma cells. Via their respective receptors, IGF-1 and IL-6 activate the RAS/ERK pathway which promotes the cell proliferation. They inhibit apoptosis through the phosphatidylinositol-3 kinase/protein kinase B/Forkhead in rabdomyosarcoma (Akt/FKHR) pathway. The cell migrates and homes to BMSCs through the SDF-1/CXCR4 axis. The IRF4 mutation, which has been associated with concomitant RAS mutations, promotes survival and proliferation. BMSCs, which are much larger cells than the myeloma cells, are shown in reduced in size in this figure as well as the receptors of both IGF-1 and IL-6. (b) The parallel evolution pattern of multiple myeloma clones resulting in intra-clonal heterogeneity. More aggressive clones result from a more aggressive N-RAS mutation in clone 2 or the acquisition of IRF4 mutation in addition to the less aggressive K-RAS in clone 4. Each clone is shown by its corresponding color in the model. (c) Snapshots of a simulation showing the competition between clones. The subclone c4 is as aggressive as the subclone c2 due to the additional IRF4 mutation and it manages to coexist with it in the remote areas with fewer cytokines. Reprinted from [Bouchnita et al., 2017a]. Color version of the picture is available on the journal website 4.1. Biology of immune response Immune processes develop in highly organized spatial structures of the lymphoid organs and the lymphatic system. We consider a part of the lymph node, i.e., the T cell zone, which contains various cell types, mainly the antigen presenting cells (APCs) and subsets of T lymphocytes. Naive T cells and some APCs (such as plasmacytoid Dendritic Cells, pDCs) enter the node with blood flow via the High Endothelial Venules (HEVs) whereas effector and/or memory T cells, and mainly DCs and macrophages home to lymph nodes via afferent lymphatic vessels [Mueller, Germain, 2009; Girard et al., 2012]. Following activation with pathogens, APCs acquire a motile state that allows their translocation to the T cell zone of draining lymph node with the afferent lymph flow [Junt et al., 2008; Förster et al., 2012]. Therefore, we assume that the influx of APCs is proportional to the level of infection in the organism. Differentiation of naive T cells into CD4+ and CD8+ T cells occurs in the thymus from progenitor T cells [Goldsby et al., 2000]. We suppose that they enter lymph nodes already differentiated and that there is a given influx of each cell type. The APCs bearing foreign antigens activate the clonal expansion of naive T lymphocytes. The activation of T cell division and death is regulated by a set of signals coming from the interactions of the antigen-specific T cell receptors (TCRs) with the MHC class I or class II presented peptides and IL-2 receptors binding IL-2. Naive T cells undergo asymmetric division [Chang, Reiner, 2008]. Some of the daughter cells continue to proliferate and differentiate. Mature CD4+ T cells produce IL-2 [Goldsby et al., 2000; Broere et al., 2011; Nelson, 2004] which influences survival and differentiation of both CD4+ and CD8+ T cells. The proliferation of CD8+ T cells is stimulated by IL-2 [Broere et al., 2011]. КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 299 Figure 6. Schematic representation of the two-scale model. Naive T cells and antigen presenting cells (APC) enter the lymph node. Due to asymmetric cell division, some T cells differentiate. Mature CD8+ T cells leave the lymph node and kill infected cells. Mature CD4+ T cells produce IL-2 that influences cell survival and differentiation. Infection level and immune cells in the body are described by ODEs. Cells in the lymph node are considered as individual objects, intracellular regulation is described by ODEs and extracellular substances by PDEs. Modified from [Bouchnita et al., 2017c] They can expand their number many thousand-fold. In addition to IL-2 enhancing the proliferation of T cells, APCs start to secrete type I Interferon (IFN) which has an antiviral- and immunomodulatory function. In fact, the effect of IFNα depends on the temporal sequence of the signals obtained by naive T cells [Welsh et al., 2012]. It can change from a normal activation of T cells followed by their proliferation and differentiation to an already differentiated state followed by apoptosis. Overall, the regulated death of T cells by apoptosis depends on the availability and the timing of TCR, IL-2 and IFN signalling. Mature CD8+ T cells (effector cells) leave the lymph node and kill infected cells. Therefore, there is a negative feedback between production of mature CD8+ T cells and the influx of APCs. In the model, an asymmetric T cell division is considered as shown in Figure 7. Naive T cell entering the draining lymph node is recruited into the immune response after the contact interaction via the T cell receptor (TCR) with APC presenting the foreign antigen. The activation and prolonged contact with APC can results in polarity of the lymphocyte. The position of the contact with the APC determines the direction of cell division and the difference between the daughter cells in terms of their differentiation state. According to [Chang, Reiner, 2008], the proximal daughter cell will undergo clonal proliferation and differentiation resulting in the generation of terminally differentiated effector cells (mature CD8+ T cells) that leave the lymph node for peripheral tissues to search and kill infected cells. The distal daughter cell becomes a memory cell. The memory cells are capable of self-renewal by slowly dividing symmetrically in the absence of recurrent infection. 4.2. Hybrid model of immune response In our model of cell dynamics, cells are considered as individual objects that can move, divide, differentiate and die. Their behavior is determined by the surrounding cells, by intracellular regulatory networks described by ordinary differential equations and by various substances in the extracellular matrix whose concentrations are described by partial differential equations. This approach was used 2019, Т. 11, № 2, С. 287–309 300 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert Figure 7. Scheme of the spatial regulation of the asymmetric T cell division in lymph nodes. Reprinted from [Bouchnita et al., 2017b] to model hematopoiesis and blood diseases [Bessonov et al., 2012]– [Stéephanou, Volpert, 2016]. Cell motion is described by equation (1.1), where ⎧ h0 − hi j ⎪ ⎪ ⎪ , h0 − hi < hi j < h0 , ⎪K ⎨ h − (h0 − h1 ) ij fi j = ⎪ ⎪ ⎪ ⎪ ⎩0, hi j ≥ h0 , hi j is the distance between the centers of the two cells i and j, h0 is the sum of their radii, K is a positive parameter and h1 is the sum of the incompressible part of each cell. The force between the particles tends to infinity if hi j decreases to h0 − h1 . Cells and concentrations We introduce the following variables. Cells in the lymph node: 1) nAPC (x, t) — the density of APCs in T cell zone; 2) nCD4 (x, t) — the density of CD4+ T cells in T cell zone (with different levels of maturity); 3) nCD8 (x, t) — the density of CD8+ T cells in T cell zone (with different levels of maturity). Extracellular variables: 4) Ie (x, t) — the concentration of IL-2 in T cell zone; 5) Ce (x, t) — the concentration of type I IFN in T cell zone. Intracellular variables: 6) Ii (t) — the intracellular concentration of IL-2-induced signalling molecules in the j-th cell; 7) Ci (t) — the intracellular concentration of type I IFN-induced signalling molecules in the j-th cell. The state variables at the level of the whole organism: 8) Ne f (t) — the total number of effector CD8+ T cells in the body; 9) Nin f (t) — the total number of infected cells in the body. КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 301 Cell division and differentiation APC and naive T cells enter the computational domain with a given frequency if there is available space. Naive T cells move in the computational domain randomly. If they contact APC, they divide asymmetrically (Figure 7). The distant daughter cell is similar to the mother cell, the proximal daughter cell becomes differentiated. When the cell reaches the half of its life cycle, it will increase its size. When it divides, two daughter cells appear, the direction of the axis connecting their centers is chosen randomly from 0 to 2π. The duration of the cell cycle is 18 hours with a random perturbation of −3 to 3 hours. We consider two levels of maturity of CD4+ T cells and three levels of CD8+ T cells. If a differentiated cell has enough IL-2 (see the next paragraph), then it divides and gives two more mature cells. Finally differentiated cells leave the lymph node. In the simulations, this means that they are removed from the computational domain. Intracellular regulation The survival and differentiation of activated CD4+ and CD8+ T lymphocytes depends on the amount of signalling via the IL-2 receptor and the type I IFN receptor. It is controlled primarily by the concentration of the above cytokines in the close proximity of the respective receptors. The signalling events lead to the up-regulation of the genes responsible for cell proliferation, differentiation and death. One can use similar type of equation to model qualitatively the accumulation of the respective intracellular signalling molecules linked to IL-2 and type I IFN receptors. The IL-2-dependent regulatory signal dynamics in individual cells can be described by the following equation: dIi α1 = Ie (xi , t) − d1 Ii . dt nT (4.1) Here Ii is the intracellular concentration of signalling molecules accumulated as a consequence of IL-2 signals transmitted through transmembrane receptor IL2R downstream the signaling pathway to control the gene expression in the j-th cell. The concentrations inside two different cells are in general different from each other. The first term in the right-hand side of this equation shows the cumulative effect of IL-2 signalling. The extracellular concentration Ie is taken at the coordinate xi of the center of the cell. The second term describes the degradation of IL-2-induced signalling molecules inside the cell. Furthermore, nT is the number of molecules internalized by T cell receptors. In a similar way, the IFN-dependent regulatory signal dynamics in individual cells can be described by the following equation: dCi α2 = Ce (xi , t) − d2 Ci . dt nT (4.2) Here Ci is the intracellular concentration of signalling molecules accumulated as a consequence of IFN signals transmitted through transmembrane receptor IFNR downstream the signaling pathway to control the gene expression in the i-th cell. The concentrations inside two different cells are in general different from each other. The first term in the right-hand side of this equation shows the cumulative effect of IFN signalling. The extracellular concentration Ce is taken at the coordinate xi of the center of the cell. The second term describes the degradation of IFN-induced signalling molecules inside the cell. To model the fate regulation of growth versus differentiation of the activated cells in relation to the timing of the IL-2 and type I IFN signalling we implement the following decision mechanism. 2019, Т. 11, № 2, С. 287–309 302 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert C1 If the concentration of activation signals induced by type I IFN, Ci , is greater than some critical level Ci∗ at the beginning of the cell cycle and that of Ii , is smaller than the critical level Ii∗ , then the cell will differentiate in a mature cell. C2 If the concentration of activation signals induced by IL-2, Ii , is greater than some critical level Ii∗ at the end of the cell cycle, then the cell will divide producing two more mature cells. C3 If Ci < Ci∗ at the beginning of cell cycle and Ii < Ii∗ at the end of cell cycle, then the cell will die by apoptosis and will be removed from the computational domain. Stochastic aspects of the model As it is discussed above, mechanical interaction of cells results in their displacement described by equation (1.1) for their centers. In order to describe random motion of cells, we add random variables to the cell velocity in the horizontal and vertical directions. Duration of cell cycle is sampled as a random variable in the interval [T − τ, T + τ]. Extracellular dynamics of cytokines. Proliferation and differentiation of T cells in the lymph node depends on the concentration of IL-2 and type I IFN. These cytokines are produced by mature CD4+ T cells and antigen-presenting cells, respectively. Their spatial distribution is described by a similar reaction-diffusion equation as follows ∂Ie = DIL ∆Ie + WIL − b1 Ie . (4.3) ∂t Here Iex is the extracellular concentration of IL-2, DIL is the diffusion coefficient, WIL is the rate of its production by CD4+ T cells, and the last term in the right-hand side of this equation describes its consumption and degradation. The production rate WIL is determined by mature CD4+ T cells. We consider each such cell as a source term with a constant production rate ρIL at the area of the cell. Let us note that we do not take into account explicitly consumption of IL-2 by immature cells in order not to introduce an additional parameter. Implicitly, this consumption is taken into account in the degradation term. For type I IFN, the equation and the terms in it have a similar interpretation: ∂Ce = DIFN ∆Ce + WIFN − b2 Ce . ∂t (4.4) Initial and boundary conditions for both concentrations IL-2 and IFN are taken zero. As before, the production rate WIFN equals ρIFN at the area filled by APC cells and zero otherwise. Infection Mature T cells leave the bone marrow. The level of CD8+ T cells (effector cells) Ne f in the body is determined by the equation dNe f = k1 T − k2 Ne f , (4.5) dt where T is their number in the lymph nodes. So the first term in the right-hand side of this equation describes production of effector cells in the lymph nodes and the second term their death in the body. Denote by Nin f the number virus-infected cells. We will describe it by the equation dNin f = f (Nin f ) − k3 Ne f Nin f . dt КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ (4.6) Hybrid models in biomedical applications 303 The first term in the right-hand side of this equation describes growth of the number of infected cells and the second term their elimination by the effector cells. The function f will be considered in the form: aNin f , f (Nin f ) = 1 + hNin f where a and h are some positive constants. Finally, the influx of APC cells into the lymph nodes is proportional to the number of infected cells Nin f . This influx is limited by the place available in the lymph node. If there is a free place sufficient to put a cell, the new cells are added. Let us also note that the lymph nodes can increase due to infection in order to produce more effector cells. 4.3. Results We illustrate the model performance by considering two scenarios, reflecting different spatial patterns of IL-2 and type IFN concentration fields. In the first one, both cytokines have the same diffusion coefficient DIL2 = DIFN , whereas in the second case the diffusion rate of IFN is 10-fold faster. The details of the numerical implementation of the hybrid model and the parameter values used for the simulations are presented in [Bouchnita et al., 2017b]. Cell population densities and cytokine concentrations are scaled with respect to some reference values. These are determined by the cell density in the lymph node ∼105 –106 mm−3 , the relative proportions of APCs , CD4+ T cells and CD8+ T cells [Ganusov, De Boer, 2007; Scandella et al., 2008; Kumar et al., 2010; Bocharov et al., 2010; Bocharov et al., 2012; Cremasco et al., 2014; Giese, Marx, 2014] and the production rate of the cytokines (described in detail in [Bouchnita et al., 2017b]). The considered cell numbers correspond to a computational domain in the T cell zone of about 100 µm × 100 µm × 100 µm. The model presented above contains two compartments, lymph node where effector cells are produced and the body where infection develops. The lymph node is described with the hybrid model while infection development in the organism by ordinary differential equations for infected cells and for effector cells. These two compartments are coupled by means of flux of effector cells from the lymph node to the body and by the flux of APC cells to the lymph node. The results of the simulations are shown in Figures 8–10. Figure 8 represents a snapshot of the lymph node T cell zone with all cells participating in the simulations: APC cells, naive T cells, differentiated CD4+ T and CD8+ T cells. Naive T cells divide when they are close to APC cells. It is an asymmetric division where a proximal daughter cell differentiates while a distant cell remains undifferentiated. Differentiated cells continue their division and maturation in the presence of IL-2 produced by mature CD4+ T cells. If the level of IL-2 is not sufficient, they die by apoptosis. Mature T cells leave the lymph node. One can see that the cytokine fields are non-uniform and their distribution patterns change essentially if the turnover parameters, e.g. the diffusion coefficient, are varied. The evolution of the number of immune cells in time depending on the ratio of diffision coefficients is shown in Figures 9 and 10 (see [Bouchnita et al., 2017c; Bouchnita et al., 2017b] for more detail). In completion, the existing studies in immunology on hybrid modelling are restricted to the consideration of immune processes on a 2D or 3D regular lattice (see Table 1 for an overview) which is a severe simplification of the physiology and anatomy of the immune system. 5. Discussion This article presents a state of the art in hybrid discrete-continuous modelling and some of its applications in biomedicine. We started by describing the underlying theory behind hybrid discretecontinuous models. In these models, cells are represented by individual objects that can move, divide, 2019, Т. 11, № 2, С. 287–309 304 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert Figure 8. Snapshot of numerical simulations of the cells and cytokines distribution in lymph node. Different cells are shown: APC (green), naive CD4+ T cells (black), naive CD8+ T cells (white), three maturity levels of differentiated CD8+ T cells (blue), two maturity levels of CD4+ T cells (yellow). Mature CD4+ T cells produce IL-2 whose concentration in the extracellular matrix is shown by the level of green. APC produce IFN (red). The upper figure shows the simulation (day 8 post infection) with equal diffusion coefficients of IL-2 and IFN, in the lower figure (day 80 post infection) the diffusion coefficient of IFN is 10 times larger than the diffusion coefficient of IL-2. Reprinted from [Bouchnita et al., 2017b]. Color version of the picture is available on the journal website Figure 9. The numbers of CD4+ and CD8+ T cells in time in the case of equal diffusion coefficients (left panel) and for the diffusion coefficient of IFN 10 times larger than the diffusion coefficient of IL-2 (right panel). Reprinted from [Bouchnita et al., 2017b] Figure 10. The numbers of APC cells (left panel) and effector T cells (right panel) in time in the case of equal diffusion coefficients (black curve) and for the diffusion coefficient of IFN 10 times larger than the diffusion coefficient of IL-2 (grey curve). Reprinted from [Bouchnita et al., 2017b] КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 305 Table 1. Hybrid models for the spatio-temporal dynamics of immune responses Hybrid model components (reference) Processes considered (duration) State space Langevin equation for T cell motility, PDE (advection-diffusion) for DCs Agent-based model for T cells [Maderazo et al., 2018] DC-based T cell vaccination 2D: deep cortical unit (100 hours) Dendritic cells (DC), T cells Agent-based model for cells, PDEs for molecules [Baldazzi et al., 2009] Clonal expansion, 3D: transport, reaction-diffusion (500 hours) DCs, B-cells, CD4+ T cells, antigen, chemokines Agent-based model for cells, ODEs for cytokines, 2D geometry of lung tissue [Fallahi-Sichani et al., 2011] Clonal expansion, 2D: chemotaxis, cell-to-cell interactions, single-cell state regulation (200 days) Macrophages, CD8+ T cells, Treg cells, Tγ cells, M. tuberculosis, TNFα, TNFR Agent-based model for cells, 3D geometry of lymph node trabecular region, location of HEVs, FRCs [Gong et al., 2013; Cilfone et al., 2015] [Marino, Kirschner, 2016] Clonal expansion, 3D: trafficking, cell-to-cell interactions (550 hrs) DCs in three states, CD4+ T cells, CD8+ T cells; Cellular Potts model for cells, PDEs for extracellular cytokines, ODEs for intracellular factors [Prokopiou et al., 2014] Clonal expansion, 3D: migration, reaction-diffusion, intracellular regulation (136 hrs) APCs, T-cells, IL-2, IL-2R, Tbet, Caspase, Fas in two states: activated/non-activated and die by apoptosis. The fate of each cell depends on the intracellular regulatory networks described by ordinary differential equations, while the extracellular regulation is modelled by partial differential equations. After the presentation of hybrid models, we have illustrated their applications in tumor growth, erythropoiesis, and immune response studies. Despite some advances in the hybrid modelling of physiological systems, there still exist important challenges that need to be overcome in order to build more reliable hybrid models. First, hybrid models are usually multiscale by construction. Hence, it can be sometimes difficult to assess the accuracy and validity of the obtained results. While it is often possible to compare the predictions of the whole model with clinical data, it is also important to validate the results of the submodels at different scales [Pathmanathan, Gray, 2018]. Secondly, hybrid models are complex which makes the interpretation of their predictions more difficult. In fact, small perturbations in the values of the input parameters sometimes leads to essential changes in the behavior of the model. Therefore, sensitivity analysis studies should be used to identify the parameters that strongly affect the behavior of the model. These studies can help in the design of more effective therapeutic strategies [Schuetz et al., 2014]. Finally, simulations with hybrid models require massive computational resources. Thus, their application is restricted to systems with a limited number of cells. However, it is possible to reduce the computational cost of a hybrid system by using efficient numerical implementation techniques such as time-stepping [Cilfone et al., 2015]. The integrative approach of hybrid models perfectly corresponds to the multiscale nature of biomedical systems. Hybrid models can be calibrated to integrate data obtained with a wide range of 2019, Т. 11, № 2, С. 287–309 306 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert acquisition techniques such as flow cytometry, nano-imaging, transcriptome sequencing analysis, etc. The increased access to biomedical data and growth of the available computational power will result in the development of more comprehensive and reliable hybrid models in the next few years. Список литературы (References) Ackleh A. S., Deng K., Ito K., Thibodeaux J. A structured erythropoiesis model with nonlinear cell maturation velocity and hormone decay rate // Mathematical biosciences. — 2006. — Vol. 204, No. 1. — P. 21–48. Anderson A., Chaplain M., Rejniak K. Single-cell-based models in biology and medicine. — Springer Science & Business Media, 2007a. Anderson A., Rejniak K., Gerlee P., Quaranta V. Modelling of cancer growth, evolution and invasion: bridging scales and models // Mathematical Modelling of Natural Phenomena. — 2007b. — Vol. 2, No. 3. — P. 1–29. Anderson A. R. A hybrid multiscale model of solid tumour growth and invasion: evolution and the microenvironment // In Single-cell-based models in biology and medicine. — Springer, 2007. — P. 3–28. Baldazzi V., Paci P., Bernaschi M., Castiglione F. Modeling lymphocyte homing and encounters in lymph nodes // BMC bioinformatics. — 2009. — Vol. 10, No. 1. — P. 387. Bauer A., Tronche F., Wessely O., Kellendonk C., Reichardt H. M., Steinlein P., Schütz G., Beug H. The glucocorticoid receptor is required for stress erythropoiesis // Genes & development. — 1999. — Vol. 13, No. 22. — P. 2996–3002. Bessonov N., Babushkina E., Golovashchenko S., Tosenberger A., Ataullakhanov F., Panteleev M., Tokarev A., Volpert V. Numerical modelling of cell distribution in blood flow // Mathematical Modelling of Natural Phenomena. — 2014. — Vol. 9, No. 6. — P. 69–84. Bessonov N., Crauste F., Fischer S., Kurbatova P., Volpert V. Application of hybrid models to blood cell production in the bone marrow // Mathematical Modelling of Natural Phenomena. — 2011. — Vol. 6, No. 7. — P. 2–12. Bessonov N., Eymard N., Kurbatova P., Volpert V. Mathematical modeling of erythropoiesis in vivo with multiple erythroblastic islands // Applied Mathematics Letters. — 2012. — Vol. 25, No. 9. — P. 1217–1221. Bessonov N., Kurbatova P., Volpert V. Dynamics of growing cell populations // Centre de Recerca Matemàtica. — 2010a. Bessonov N., Kurbatova P., Volpert V. Particle dynamics modelling of cell populations // Mathematical Modelling of Natural Phenomena. — 2010b. — Vol. 5, No. 7. — P. 42–47. Bessonov N., Pujo-Menjouet L., Volpert V. Cell modelling of hematopoiesis // Mathematical Modelling of Natural Phenomena. — 2006. — Vol. 1, No. 2. — P. 81–103. Bessonov N., Reinberg N., Volpert V. How morphology of artificial organisms influences their evolution // Ecological complexity. — 2015. — Vol. 24. — P. 57–68. Bocharov G., Chereshnev V., Gainova I., Bazhan S., Bachmetyev B., Argilaguet J., Martinez J., Meyerhans A. Human immunodeficiency virus infection: from biological observations to mechanistic mathematical modelling // Mathematical Modelling of Natural Phenomena. — 2012. — Vol. 7, No. 5. — P. 78–104. Bocharov G., Züst R., Cervantes-Barragan L., Luzyanina T., Chiglintsev E., Chereshnev V. A., Thiel V., Ludewig B. A systems immunology approach to plasmacytoid dendritic cell function in cytopathic virus infections // PLoS pathogens. — 2010. — Vol. 6, No. 7. — P. e1001017. КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 307 Bouchnita A., Belmaati F.-E., Aboulaich R., Koury M. J., Volpert V. A hybrid computation model to describe the progression of multiple myeloma and its intra-clonal heterogeneity // Computation. — 2017a. — Vol. 5, No. 1. — P. 16. Bouchnita A., Bocharov G., Meyerhans A., Volpert V. Hybrid approach to model the spatial regulation of T cell responses // BMC immunology. — 2017b. — Vol. 18, No. 1. — P. 29. Bouchnita A., Bocharov G., Meyerhans A., Volpert V. Towards a multiscale model of acute hiv infection // Computation. — 2017c. — Vol. 5, No. 1. — P. 6. Bouchnita A., Eymard N., Moyo T. K., Koury M. J., Volpert V. Bone marrow infiltration by multiple myeloma causes anemia by reversible disruption of erythropoiesis // American journal of hematology. — 2016. — Vol. 91, No. 4. — P. 371–378. Brioli A., Melchor L., Cavo M., Morgan G. J. The impact of intra-clonal heterogeneity on the treatment of multiple myeloma // British journal of haematology. — 2014. — Vol. 165, No. 4. — P. 441–454. Broere F., Apasov S. G., Sitkovsky M. V., van Eden W. A2 T cell subsets and T cell-mediated immunity // In Principles of immunopharmacology. — Springer. — P. 15–27. Chang J., Reiner S. Asymmetric division and stem cell renewal without a permanent niche: lessons from lymphocytes // In Cold Spring Harbor symposia on quantitative biology. — 2008. — Vol. 73. — P. 73–79. — Cold Spring Harbor Laboratory Press. Chasis J. A., Mohandas N. Erythroblastic islands: niches for erythropoiesis // Blood. — 2008. — Vol. 112, No. 3. — P. 470–478. Cilfone N. A., Kirschner D. E., Linderman J. J. Strategies for efficient numerical implementation of hybrid multi-scale agent-based models to describe biological systems // Cellular and molecular bioengineering. — 2015. — Vol. 8, No. 1. — P. 119–136. Crauste F., Demin I., Gandrillon O., Volpert V. Mathematical study of feedback control roles and relevance in stress erythropoiesis // Journal of theoretical biology. — 2010. — Vol. 263, No. 3. — P. 303–316. Crauste F., Pujo-Menjouet L., Génieys S., Molina C., Gandrillon O. Adding selfrenewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis // Journal of theoretical biology. — 2008. — Vol. 250, No. 2. — P. 322–338. Cremasco V., Woodruff M. C., Onder L., Cupovic J., Nieves-Bonilla J. M., Schildberg F. A., Chang J., Cremasco F., Harvey C. J., Wucherpfennig K. et al. B cell homeostasis and follicle confines are governed by fibroblastic reticular cells // Nature immunology. — 2014. — Vol. 15, No. 10. — P. 973. De Maria R., Testa U., Luchetti L., Zeuner A., Stassi G., Pelosi E., Riccioni R., Felli N., Samoggia P., Peschle C. Apoptotic role of fas/fas ligand system in the regulation of erythropoiesis // Blood. — 1999. — Vol. 93, No. 3. — P. 796–803. Drasdo D. Center-based single-cell models: An approach to multi-cellular organization based on a conceptual analogy to colloidal particles // In Single-Cell-Based Models in Biology and Medicine. — Springer, 2007. — P. 171–196. Fallahi-Sichani M., El-Kebir M., Marino S., Kirschner D. E., Linderman J. J. Multiscale computational modeling reveals a critical role for tnf-α receptor 1 dynamics in tuberculosis granuloma formation // The Journal of Immunology. — 2011. — P. 1003299. Fischer S., Kurbatova P., Bessonov N., Gandrillon O., Volpert V., Crauste F. Modeling erythroblastic islands: using a hybrid model to assess the function of central macrophage // Journal of theoretical biology. — 2012. — Vol. 298. — P. 92–106. Förster R., Braun A., Worbs T. Lymph node homing of T cells and dendritic cells via afferent lymphatics // Trends in immunology. — 2012. — Vol. 33, No. 6. — P. 271–280. Gandrillon O. The v-erba oncogene // In Thyroid Hormone Receptors. — Springer, 2002. — P. 91–107. 2019, Т. 11, № 2, С. 287–309 308 N. M. Bessonov, G. A. Bocharov, A. Bouchnita, V. A. Volpert Gandrillon O., Schmidt U., Beug H., Samarut J. TGF-β cooperates with TGF-α to induce the self– renewal of normal erythrocytic progenitors: evidence for an autocrine mechanism // The EMBO journal. — 1999. — Vol. 18, No. 10. — P. 2764–2781. Ganusov V. V., De Boer R. J. Do most lymphocytes in humans really reside in the gut? // Trends in immunology. — 2007. — Vol. 28, No. 12. — P. 514–518. Giese C., Marx U. Human immunity in vitro–solving immunogenicity and more // Advanced drug delivery reviews. — 2014. — Vol. 69. — P. 103–122. Girard J.-P., Moussion C., Förster R. HEVs, lymphatics and homeostatic immune cell trafficking in lymph nodes // Nature Reviews Immunology. — 2012. — Vol. 12, No. 11. — P. 762. Goldsby R., Kuby J., Kindt T. Immunology. — WH Freeman & Co (Sd), 2000. Gong C., Mattila J. T., Miller M., Flynn J. L., Linderman J. J., Kirschner D. Predicting lymph node output efficiency using systems biology // Journal of theoretical biology. — 2013. — Vol. 335. — P. 169–184. Junt T., Scandella E., Ludewig B. Form follows function: lymphoid tissue microarchitecture in antimicrobial immune defence // Nature Reviews Immunology. — 2008. — Vol. 8, No. 10. — P. 764. Karttunen M., Vattulainen I., Lukkarinen A. Novel methods in soft matter simulations. — Vol. 640. — Springer Science & Business Media, 2004. Koulnis M., Liu Y., Hallstrom K., Socolovsky M. Negative autoregulation by fas stabilizes adult erythropoiesis and accelerates its stress response // PLoS One. — 2011. — Vol. 6, No. 7. — P. e21192. Koury M. J., Bondurant M. C. Erythropoietin retards dna breakdown and prevents programmed death in erythroid progenitor cells // Science. — 1990. — Vol. 248 (4953). — P. 378–381. Kumar V., Scandella E., Danuser R., Onder L., Nitschké M., Fukui Y., Halin C., Ludewig B., Stein J. V. Global lymphoid tissue remodeling during a viral infection is orchestrated by a b cell-lymphotoxindependent pathway // Blood. — 2010. — Vol. 115. — P. 4725–4733. Kurbatova P., Bernard S., Bessonov N., Crauste F., Demin I., Dumontet C., Fischer S., Volpert V. Hybrid model of erythropoiesis and leukemia treatment with cytosine arabinoside // SIAM Journal on Applied Mathematics. — 2011. — Vol. 71, No. 6. — P. 2246–2268. Maderazo D. L., Flegg J. A., Neeland M. R., de Veer M. J., Flegg M. B. Physiological factors leading to a successful vaccination: A computational approach // Journal of theoretical biology. — 2018. — Vol. 454. — P. 215–230. Mahaffy J. M., Bélair J., Mackey M. C. Hematopoietic model with moving boundary condition and state dependent // J. theor. Biol. — 1998. — Vol. 190. — P. 135–146. Marino S., Kirschner D. E. A multi-compartment hybrid computational model predicts key roles for dendritic cells in tuberculosis infection // Computation. — 2016. — Vol. 4, No. 4. — P. 39. Melchor L., Brioli A., Wardell C., Murison A., Potter N., Kaiser M., Fryer R., Johnson D., Begum D., Wilson S. H. et al. Single-cell genetic analysis reveals the composition of initiating clones and phylogenetic patterns of branching and parallel evolution in myeloma // Leukemia. — 2014. — Vol. 28, No. 8. — P. 1705. Morgan G. J., Walker B. A., Davies F. E. The genetic architecture of multiple myeloma // Nature Reviews Cancer. — 2012. — Vol. 12, No. 5. — P. 335. Mueller S. N., Germain R. N. Stromal cell contributions to the homeostasis and functionality of the immune system // Nature Reviews Immunology. — 2009. — Vol. 9, No. 9. — P. 618. Nelson B. H. Il-2, regulatory T cells, and tolerance // The Journal of Immunology. — 2004. — Vol. 172, No. 7. — P. 3983–3988. КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Hybrid models in biomedical applications 309 Pain B., Woods C., Saez J., Flickinger T., Raines M., Peyroll S., Moscovici C., Moscovici M., Kung H.-J., Jurdic P. et al. EGF-r as a hemopoietic growth factor receptor: the c-erbb product is present in chicken erythrocytic progenitors and controls their self-renewal // Cell. — 1991. — Vol. 65, No. 1. — P. 37–46. Palumbo A., Anderson K. Multiple myeloma // N. Engl. J. Med. — 2011. — Vol. 364, No. 1. — P. 51–61. Pathmanathan P., Gray R. A. Validation and trustworthiness of multiscale models of cardiac electrophysiology // Frontiers in Physiology. — 2018. — Vol. 9. — P. 106. Prokopiou S. A., Barbarroux L., Bernard S., Mafille J., Leverrier Y., Arpin C., Marvel J., Gandrillon O., Crauste F. Multiscale modeling of the early CD8 T-cell immune response in lymph nodes: an integrative study // Computation. — 2014. — Vol. 2, No. 4. — P. 159–181. Rubiolo C., Piazzolla D., Meissl K., Beug H., Huber J. C., Kolbus A., Baccarini M. A balance between raf-1 and fas expression sets the pace of erythroid differentiation // Blood. — 2006. — Vol. 108, No. 1. — P. 152–159. Sakr S., Jeanjean R., Zhang C.-C., Arcondeguy T. Inhibition of cell division suppresses heterocyst development in anabaena sp. strain pcc 7120 // Journal of bacteriology. — 2006. — Vol. 188, No. 4. — P. 1396–1404. Sawyer S. T., Jacobs-Helber S. M. State-of-the-art review: Unraveling distinct intracellular signals that promote survival and proliferation: Study of erythropoietin, stem cell factor, and constitutive signaling in leukemic cells // Journal of hematotherapy & stem cell research. — 2000. — Vol. 9, No. 1. — P. 21–29. Scandella E., Bolinger B., Lattmann E., Miller S., Favre S., Littman D. R., Finke D., Luther S. A., Junt T., Ludewig B. Restoration of lymphoid organ integrity through the interaction of lymphoid tissue–inducer cells with stroma of the t cell zone // Nature immunology. — 2008. — Vol. 9, No. 6. — P. 667. Schuetz T. A., Mang A., Becker S., Toma A., Buzug T. M. Identification of crucial parameters in a mathematical multiscale model of glioblastoma growth // Computational and Mathematical Methods in Medicine. — 2014. Spivak J. L., Pham T., Isaacs M., Hankins W. D. Erythropoietin is both a mitogen and a survival factor // Blood. — 1991. — Vol. 77, No. 6. — P. 1228–1233. Stéephanou A., Volpert V. Hybrid modelling in biology: a classification review // Mathematical Modelling of Natural Phenomena. — 2016. — Vol. 11, No. 1. — P. 37–48. Tosenberger A., Ataullakhanov F., Bessonov N., Panteleev M., Tokarev A., Volpert V. Modelling of platelet–fibrin clot formation in flow with a dpd–pde method // Journal of mathematical biology. — 2016. — Vol. 72, No. 3. — P. 649–681. Tsiftsoglou A. S., Vizirianakis I. S., Strouboulis J. Erythropoiesis: model systems, molecular regulators, and developmental programs // IUBMB life. — 2009. — Vol. 61, No. 8. — P. 800–830. Welsh R. M., Bahl K., Marshall H. D., Urban S. L. Type 1 interferons and antiviral CD8 T-cell responses // PLoS pathogens. — 2012. — Vol. 8, No. 1. — P. e1002352. Wichmann H., Loeffler M., Pantel K., Wulff H. A mathematical model of erythropoiesis in mice and rats part 2: Stimulated erythropoiesis // Cell Proliferation. — 1989. — Vol. 22, No. 1. — P. 31–49. 2019, Т. 11, № 2, С. 287–309