# Mathematical model of T-cell lymphoblastic lymphoma: disease, treatment, cure or relapse of a virtual cohort of patients

Mathematical model of T-cell lymphoblastic lymphoma: disease, treatment, cure or relapse of a... Abstract T lymphoblastic lymphoma (T-LBL) is a rare type of lymphoma with a good prognosis with a remission rate of 85%. Patients can be completely cured or can relapse during or after a 2-year treatment. Relapses usually occur early after the remission of the acute phase. The median time of relapse is equal to 1 year, after the occurrence of complete remission (range 0.2–5.9 years) (Uyttebroeck et al., 2008). It can be assumed that patients may be treated longer than necessary with undue toxicity. The aim of our model was to investigate whether the duration of the maintenance therapy could be reduced without increasing the risk of relapses and to determine the minimum treatment duration that could be tested in a future clinical trial. We developed a mathematical model of virtual patients with T-LBL in order to obtain a proportion of virtual relapses close to the one observed in the real population of patients from the EuroLB database. Our simulations reproduced a 2-year follow-up required to study the onset of the disease, the treatment of the acute phase and the maintenance treatment phase. 1. Introduction This work was carried out as a part of the CRESim study (Child-Rare-Euro-Simulation), a research project funded by the ERA-NET PrioMedChild (Priority Medicines for Children). The main objective of the CRESim project was to create a platform for performing in silico experiments assessing and optimizing designs of RCTs for drug evaluation in children with rare diseases. Virtual RCTs can be simulated using different experimental designs and compared in terms of trial duration (for the patient, the investigator and the sponsor) and precision of the estimation of the treatment effect. In this particular work, the CRESim methodology was applied to study the duration of the maintenance treatment of T-cell lymphoblastic lymphoma. Lymphoblastic lymphoma (LBL) is a rare form of neoplasm that can develop from immature precursor cells in the bone marrow (BM) or from thymic T cells at varying stages of differentiation (Cortelazzo et al. (2011); Schmidt & Burkhardt (2013)). In developed countries, the annual cancer incidence in children under 15 years of age is estimated at 140 new cases per million and the LBL annual incidence at 0.3 to 0.5 cases per 100,000 children and adolescents (Brugiéres et al. (2012); Hochberg & Cairo (2009)). T lymphoblastic lymphoma (T-LBL) is a rare disease with few mathematical models developed for this type of cancer. A competition between normal and tumourous clonotypes of naive T-cells has been modelled in Chrobak et al. (2012) as a continuous-time bivariate Markov process. The study presented in Page & Uhr (2005) illustrates how the interaction between tumour and antibody for the murine BCL1 lymphoma leads to the dormancy of the tumour. The model takes into account cell cycle arrest and apoptosis of the tumour cells due to the immune response. The model of Non-Hodgkins lymphomas (NHL) (Ribba et al., 2005) includes a tumour structure incorporating mature and immature vessels, vascular structural adaptation and NHL cell-cycle kinetics in addition to pharmacokinetics and pharmacodynamics of the used drug (Doxorubicin). Another model of non-Hodgkin’ s lymphoma considers the tumour as a mixture of cells, interstitial fluid, and extracellular matrix (Frieboes et al., 2013). The concentrations of cells, as well as of oxygen and nutrients are described by partial differential equations. A recent review of lymphoma and leukaemia models (Clapp & Levy, 2014) focuses on haematopoiesis and drug improvement studies, and especially on drug resistance. Most T-LBL patients typically present with a mediastinal tumour. We therefore restrict our focus on a model of normal thymus, modelling three phases: its invasion by a tumour due to proliferation of a single cancer cell, the treatment of the acute phase and the maintenance treatment phase. For validation purposes, the results of our simulations compare the proportions of relapses in the virtual population with the ones from the EuroLB database of real patients. The aim is to build a model allowing the variation of the proportion of relapses of a population of T-LBL patients as a function of the drug dose and as a function of the treatment duration. The simulations should allow to test each of these essential elements separately or together. Inter-patients variability needed to reproduce different behaviours, cured or not, of real patients can be simulated by adding random noise to several existing models. However, to reach our objective, a more complex model is required. We use a mathematical model at a cellular level to simulate tumours and healthy thymuses; the presence or not, of surviving cancer cells after chemotherapy in a virtual patient determines his state of health. Given the clinical presentation of the disease, our model assumes a tumour arising in the thymus. We begin with the simulation of a healthy thymus. The model included the onset of the disease, the treatment of its acute and maintenance phase. In accordance with biomedical data, we assume that patient response to treatment depends on sensitivity of cancer cell to the intracellular drug concentration. In the first part of the work, we model the proliferation of a malignant tumour at a cell-scale. In this model, a given sensitivity of cancer cells to treatment is assigned to each cell as a parameter. Then, taking into account the probability density function of treatment sensitivity, the therapeutic effect of drug is simulated for virtual patients. 2. Thymus modelling The thymus is an important organ in the development of the immune system in young children. Its function becomes negligible in adults. The thymus is composed of several lobules, containing a cortex and a medulla separated by a cortico medullary junction. The thymus, in particular the cortex, is the site of maturation Brousse (2009), education and selection of T lymphocytes. Our simulations are focused on the processes called ontogenesis. 2.1 Healthy functioning thymus modelling: Ontogenesis The mathematical model of T cells ontogenesis takes into account four types of T cells: precursors, double negative (DN), double positive (DP) and simple positive (SP). The mechanisms allowing their progression between the consecutive steps implemented in our model are summarized in Fig. 1. The given percentages of apoptosis and the proportion of each kind of T cells are observed in reality and are computed by the model. In a healthy thymus, the proportion of T cells remains approximately constant: 3% of DN, 85% of DP and 12% of SP cells. Fig. 1. View largeDownload slide Discursive representation of the simulated ontogenesis. Fig. 1. View largeDownload slide Discursive representation of the simulated ontogenesis. Each day, $$5 \times 10^7$$ mature cells are produced but only $$10^6$$ cells leave the thymus. 98% of the thymocytes die by apoptosis by $$\beta$$, positive and negative selections. During their development, lymphocytes express several receptors allowing them to receive signals from their environment. Cells that fail to produce a functional pre-TCR (T Cell Membrane Receptor) are eliminated by apoptosis. This is called the $$\beta$$-selection. In positive selection, through contact with cortical epithelial cells, DN lymphocytes receive a signal of life or death. In negative selection, through contact with epithelial cells of medulla, DP surviving lymphocytes also receive this signal. The selection processes depend on the affinity TCRs for self MHC (Major Histocompatibility Complex) molecules from epithelial cells. Apoptosis occurs when there is a low affinity for positive selection and a high affinity for negative selection. Stem cells in the BM differentiate into thymic precursors and enter the thymus through blood vessels. After several stages of differentiation and proliferation, they then become T cells. Precursors are committed to a specific lineage through signalling. We only consider the T lineage by activation of their receptors $$Notch$$. During the early stages of maturation, precursors become DN1s, which are multipotent cells able to differentiate into T lymphocytes, NK (Natural Killer) lymphocytes and dendritic cells. Pro T thymocytes (DN2) mature and become DPs or NKs. Finally, DPs become SPs and exit the thymus through blood vessels (Janeway et al., 2011). This classification depends on the expression of their receptors CD4 and CD8, DP: $$CD4^+CD8^+$$, DN: $$CD4^-CD8^-$$, SP: $$CD4^+CD8^-$$ or $$CD4^-CD8^+$$. The maturation, migration, proliferation and differentiation of T cells lineage are controlled by interactions with other cells in the thymus, especially with cortical and medullary epithelial cells, nurse cells and dendritic cells. Our model only takes into account their actions on DNs, DPs and SPs and not their maturation processes. In the normal thymus, mature cells produce cytokines such as interleukins, growth factors and inhibitory factors. They provide a biological support for proliferation, differentiation or death signal. It is suggested that intracellular signals through different MAPK (mitogen-activated protein kinase) cascades selectively guide positive and negative selection of T lymphocytes (Hoelzer & Gokbuget, 2009). They are regulators of the differentiation of immature thymocytes from DN to DP cells, most probably acting as a transducer for pre-T cell receptor signalling (Crompton et al., 1996). The ERK (Extracellular Signal-regulated Kinease) and the p38 (a Mitogen-activated protein kinease) signalling cascades are also critically involved in TCR signals inducing positive selection and negative selection, respectively (Crompton et al., 1996). Differentiation, proliferation and self-renewal of real T cells are controlled by signals emerging from their microenvironment. In our model, the balance between cytokines families with antagonist roles is summarized by a reaction diffusion equation and by a death probability. At each differentiation stage a high rate of apoptosis is observed mostly due to positive and negative selections. Apoptosis of T cells, whatever the cause, is modelled by a probability of remaining alive for each cell type (cf. Subsection 3.1.4). Proliferation and differentiation of thymocytes are responses to signals sent by epithelial cells modelled by an extra regulation using partial differential equation (Subsections 3.1.2) that influences the T cells intraregulation simulated by an ordinary differential equation (cf. Subsection 3.1.3). 2.2 Disease modelling The disease begins in the thymus but the tumour and the effect of the treatment can also reach the whole body through blood vessels and affect the BM. The tumour results from a malignant transformation of clone abnormal progenitor cells with infinite self-renewal potential which occurs in the thymus (MA et al., 2002). Such cells present a genetic lesion that causes their proliferation and blocks their differentiation (Chiaretti & Foà R., 2009). Genetic aberrations in leukaemogenesis can affect transcription factor, protein-protein interaction, signal transduction, fate determination and cell cycle activator and differentiation (Chiaretti & Foà R., 2009). The consequences of these genetic mutations are very complex and are not modelled in our study. In 90% of cases, lymphomas are caused by an excessive proliferation of immature thymocytes Brousse (2009). Cancer cells appear to proliferate rather than to differentiate, causing the emergence and development of a tumour. Consequently, the tumour is modelled by a set of cells that proliferate with a low rate of apoptosis (probability of death close to zero) and without any differentiation. The intra and extraregulations, which induce differentiation, are not taken into account for the simulation of cancer cells. 3. Mathematical and computational model We visualize the simplified operating process of the thymus during two years using an hybrid model with an interface made with $$C^{++}$$ graphical library wxWidget (Fischer et al., 2012) included in specific software developed for haematopoiesis modelling (Bessonov et al., 2009). In this computed 2D representation, T cells are represented as circles with specific colours for each stage of differentiation, (Figs 2 and 3). Cell divisions are model as follow: during cell cycle, cell radius grows linearly. It becomes twice the initial radius at the end of the cycle. Then, each mother cell is replaced by two daughter cells. Geometrically, two circles with the initial radius are located inside the circle with a twice larger radius. The direction of cell division, that is the direction of the line connecting the centres of daughter cells is random (Fischer et al., 2012). Fig. 2. View largeDownload slide Development of the tumour in the thymus. The number of malignant cells grows exponentially with time. They invade the thymus and push out normal cells. The volume of the thymus increases and gradually cancer cells replace the normal T-cells. Cells shown in the figure are: DN (blue), DP (green), SP (red), stromal cells (large purple cells), malignant cells M (purple). The extracellular substances produced by nurse cells are shown as green halos. Simulations made with the basic parameters values. (Color code refers to the electronic version of the paper.) Fig. 2. View largeDownload slide Development of the tumour in the thymus. The number of malignant cells grows exponentially with time. They invade the thymus and push out normal cells. The volume of the thymus increases and gradually cancer cells replace the normal T-cells. Cells shown in the figure are: DN (blue), DP (green), SP (red), stromal cells (large purple cells), malignant cells M (purple). The extracellular substances produced by nurse cells are shown as green halos. Simulations made with the basic parameters values. (Color code refers to the electronic version of the paper.) Fig. 3. View largeDownload slide Simulations of lymphoma treatment by chemotherapy. Cancer cells are killed by the treatment. When the tumour disappears, precursors enter the thymus and recreate its normal structure after several stages of differentiation and proliferation. Precursors (gray), DN (blue), DP (green), SP (red), stromal cells (large purple cells) malignant cells (purple). Simulations made with basic parameters values with the exception of drug dose of acute treatment, equal to 5 (see explanation of this choice in the text). (Color code refers to the electronic version of the paper.) Fig. 3. View largeDownload slide Simulations of lymphoma treatment by chemotherapy. Cancer cells are killed by the treatment. When the tumour disappears, precursors enter the thymus and recreate its normal structure after several stages of differentiation and proliferation. Precursors (gray), DN (blue), DP (green), SP (red), stromal cells (large purple cells) malignant cells (purple). Simulations made with basic parameters values with the exception of drug dose of acute treatment, equal to 5 (see explanation of this choice in the text). (Color code refers to the electronic version of the paper.) On a computer view, a cell is a set of parameters (the meaning of all the considered parameters and their values are given in Appendix A.2). In order not to have a predetermine cell behaviour, these values can vary with time and cells location. Each cell has its own initial and critical values and coefficients of equations considered in the model. Their fate depend on the computation results done during their whole life (cf. Section 3.2). To increase legibility of this representation, cells cannot overlap. It is therefore necessary to consider their positions and displacements (Subsection 3.1.1, equations (3.4) and (3.3). A more complete approach of cells motion in hybrid models can be found in Kurbatova et al. (2013); Eymard et al. (2015); Eymard & Kurbatova (2015); Bessonov et al. (2011); Fischer et al. (2012). Displacements of cells in the thymus are important as it corresponds to a biologic reality; their maturation depends on their location in the thymus. The equations of the model are solved with an algorithm involving the spatial discretization of a computational domain, that is the considered thymus cut, and discretization of the time of simulation, that is each iteration of the computation corresponds to a time step $$\Delta t$$. In numerical simulations, equation (3.6) is solved by a finite difference method with Thomas algorithm and an alternative direction method. Neumann (no-flux) boundary conditions are considered as the boundary of the rectangular domain (Eymard et al., 2015; Fischer et al., 2012). Constant space and time steps are used. Equations (3.4) and (3.5) are solved by Euler method. The outputs of our model are the amounts of each kind of cells at any time. The results obtained with this approach cannot be directly compared with medical data. However, to our knowledge no medical examinations are able to provide such data. We defined the number of relapses at the end of the therapeutic maintenance phase simulation as follow: if the number of cancer cells is equal to zero, the patient is cured, otherwise the patient relapses. 3.1 T cell ontogenesis simulation 3.1.1 Hybrid discrete-continuous models Hybrid discrete-continuous models are widely used in the investigation of dynamics of cell populations in biological tissues and of organisms, which involve processes at different scales. This model is used to describe the education of T cells, the growth of the tumour and the action of the treatment on tumour cells. Hybrid models consist in coupling two models: one with a discrete approach for healthy and cancer cells modelled as objects interacting with their environment, and another with a continuous approach using ordinary and partial differential equations to describe intracellular and extracellular regulations that determine the fate of each cell (Alarcón et al., 2003; Anderson et al., 2007; Anderson & Chaplain, 1998; Dallon, 2007; Deutsch, 2007; Eymard et al., 2015; Kurbatova et al., 2013; Fischer et al., 2012). We consider the following equation to describe intracellular proteins in one cell:   dui(t)dt=F(ui(t),v(xi,t)), (3.1) where $$u_i$$ is a vector of intracellular concentrations of the cell $$i$$, $$v$$ is a vector of extracellular concentrations, $$F$$ is the vector of reaction rates which will be specified below and $$x_i$$ is the vector of coordinates of the centre of the cell. The quantity of protein depends on the concentrations of nutrients or biological products in the extracellular matrix described by the reaction-diffusion equation:   ∂v∂t=DΔv+H(v,c), (3.2) where $$c$$ is the local cell density, $$D$$ is the diffusion coefficient and $$H$$ is the rate of consumption or production of these substances by cells. These compounds can either be nutrients coming from outside or some other bio-chemical products consumed or produced by cells (Kurbatova et al., 2013; Eymard et al., 2015; Eymard & Kurbatova, 2015; Bessonov et al., 2011; Fischer et al., 2012; Kurbatova et al., 2011). The value of extracellular concentration $$v(x_i,t)$$ in equation (3.1) is taken at $$x=x_i$$. In our model, when cells divide, they push each other and cause their displacement inside the thymus. In order to describe mechanical interactions between cells, we restrict ourselves to the simplest model where cells are represented as compressible elastic spheres with an incompressible core (black circle inside each cell, Figs 2 and 3). Cell motion is determined by the sum of mechanical forces acting on a cell by other cells. Under the assumption of small deformations, the mechanical contact forces acting between cells can be expressed as a function $$f(d_{ij})$$ of the distance $$d_{ij}$$ between their centres. Cell motion is described by the displacement of its centre by Newton’s second law (Eymard & Kurbatova, 2015; Fischer et al., 2012):   mxi¨+μmxi˙−∑j≠if(dij)=0, (3.3) where $$m$$ is the mass of the cell, $$x_i$$ is the centre of the $$i$$-th cell, the second term in the left-hand side of this equation describes the friction by the surrounding medium, the last term represents the sum of potential forces acting on the i-th cell from all other cells (Bessonov et al., 2011; Fischer et al., 2012); (Kurbatova et al., 2011). The force between two spherical particles is given by the equation:   f(dij)={+∞,dij<d0−2H1κd0−dijdij−d0+2H1,dij<d0,0,dij≥d0  (3.4) where $$d_0$$ is the sum of cells radii, $$\kappa$$ is a positive parameter, and $$H_1$$ is the size of the outer compressible part of the cell. The force between the particles tends to infinity when $$d_{ij}$$ decreases to $$d_0-2H_1$$. The force is repulsive if $$d_{ij} < d_0$$ and is null otherwise. 3.1.2 Intracellular regulation of T cell In vivo experiments show that growth factors promote proliferation of early thymocyte cells (Martins et al., 2012; Peaudecerf et al., 2012). During their maturation, thymocyte cells move and away from nurse cells and they continue to differentiate. In order to describe the self-renewal and differentiation of thymocytes, we introduce an intracellular variable $$u$$ which corresponds to the concentration of the proteins MAPK. Its evolution is described by the equation:   dudt =c0+kG, (3.5) where $$c_0$$ is a positive constant and $$G$$ is the concentration of extracellular growth factors. If $$u$$ exceeds some critical value $$u_ {crit}$$ at the end of cell cycle, then the cell self-renews, otherwise it differentiates. 3.1.3 Extracellular regulation of T cell Early T cells are first located in a stromal cell niche and exposed to stromal factors (Boehm, 2012), which determine their proliferation and differentiation (Crompton et al., 1996; Rincon et al., 2001). Nurse cells produce growth factor $$G$$ which promotes self-renewal of thymocytes. If thymocyte cells are close to nurse cells, their intracellular concentration, denoted $$u$$, will increase due to the presence of growth factor. It will reach the critical value and the cell will self-renew. The cells located sufficiently far away from the nurse cells will differentiate. The concentration of the growth factor $$G$$ is described by the reaction-diffusion equation:   ∂G∂t=DΔG+W−γG, (3.6) where $$D$$ is the diffusion coefficient, $$W$$ is the production rate and $$\gamma$$ is the degradation rate blue of growth factor. These parameters are assumed to be constant. The values of parameters and initial conditions described in Subsections 3.1.2 and 3.1.3 are specified in Table 2 of Appendix A (Subsection A.2). This Table gives the basic value of parameters used for all the presented simulations unless otherwise specified in their legends. 3.1.4 $$\beta$$-selection, negative and positive selections Signal transductions are responsible for the different selections of T cells. Epithelial cells form a fibrous network in the thymus (Brousse, 2009). We assume that they are numerous enough and uniformly distributed in space such that each T cell can receive their signals of life and death from a cell to cell contact. In our simulations, apoptosis is modelled by a probability of remaining alive for each cell type. The probability of survival is equal to 60% for precursors, 40% for DN, 50% for DP and 40% for SP. After all maturational stages, the probability of 95% of apoptosis for healthy simulated cells is observed, this ratio is close to the ratio obtained by biological cells. 3.2 Thymus simulation The simulated cut of thymus is a set of cells, characterized by a list of parameters values. Initial values are assigned to cells at birth, such as: life expectancy $$T$$ (the term ‘cell cycle’ will also be used to define this duration), and initial concentration of protein or critical values such as lethal quantity of drug $$p_{crit}$$ and needed amount of intracellular protein for self renewal or differentiation $$u_{crit}$$. These parameters can be transmitted from mother to daughter cells or they can be random variables. At each iteration of the computation and for each cell, the position, the intracellular amount of drug $$p$$, of proteins $$u$$ and of growth factor available are updated, by solving equations (3.3, 3.4, 4.2, 3.5 and 3.6), respectively. It is important to note that most of the parameters are not biological characteristics of cells. They summarize the causes of a behaviour. Several proteins, growing factors and MAPK are involved in the choice of cell between apoptosis, self renewal and differentiation. In our model, a single protein $$u$$ and a single extracellular substance $$G$$ are considered. Similarly, the values of $$u_{crit}$$, $$k$$ and $$u_0$$ are chosen to obtain division rates close to the reality. The choice of the parameters’ values is a result of the calibration of the model by retrospective comparison of computational outputs with known real data. We compared the number of relapses after a two-year treatment observed in the real EuroLB patients and the virtual patients of our model. All other factors being equal, we varied the duration of treatments and the doses of therapy (cf. Section 4, results shown in Fig. 6). This can be considered as a sensitivity analysis of the most important inputs of the model: the duration of the treatment and its dose. However, a sensitivity analysis has not been performed for all others parameters, these simplifications are necessary to keep a reasonable computation time (the simulation of two patients requires about one day of the CPU time -bi-processeurs Intel Xeon 5650- of 10 nodes of computer cluster, namely more than 3 months of computations for all simulations). Cell cycle duration is taken in average of $$T=3$$ days plus/minus a random variation $$\delta$$. We set $$\delta = 1$$. Cell cycle varies in the time interval $$[T-\delta, T+\delta]$$, with a uniform distribution. After a period of time equal to its life expectancy, the cell chooses between self-renewal, differentiation or apoptosis. In the model, the amounts of intracellular proteins for healthy cells and of drug for cancer cells are compared with the corresponding threshold values and a decision is taken between the three options. When cells become SPs, they leave the thymus. To maintain a constant number of each cell types, new cells are introduced in the computational domain if there is some space available for them, by a periodic input of precursors. Consequently, new cells appear if the number of lymphocytes in the thymus is sufficiently low. In the case of tumour growth, the number of cancer cells is high and no new cells appear in the thymus. Spatial cell organization simulations with different cell types are shown in Fig. 2.Figure A.1 shows the results of numerical simulations of healthy thymus which correspond to one year of real time. The number of the different cell types is constant in average with some oscillations related to cell cycle. 4. Lymphoma development and treatment 4.1 Simulation of individual behaviour of a patient with a lymphoma: the onset of the disease and its treatment. 4.1.1 Simulation of spontaneous evolution of lymphoma The origin of a tumour is related to early thymocytes which have an excessive self-renewal potential and insufficient differentiation and apoptosis because of genetic mutations. To mimic tumour onset, we introduce a single mutated cell which proliferates. We model a tumour as a set of cancer cells. Malignant cells self-renew at the end of each cell cycle. They have a very low value of $$u_{crit}$$ and stay alive with a given probability $$q$$ ($$q=99\%$$ in all the simulations). In order to model the action of chemotherapy treatment on cancer cells, we prescribe to each cell the value $$p_{crit}$$ of the intracellular drug concentration needed to kill this particular cell. We characterize the cell population by an average value $$p_{crit}^m$$ with a normal distribution around this average value with standard deviation $$\sigma$$ equal 10% of the average value. In the computational domain, the multiplication of cancer cells increases the volume of the thymus and gradually invades it entirely. Healthy cells cannot develop normally and are progressively replaced by malignant cells. The tumour invades the thymus and destroys the spatial cell organization which is necessary for normal maturation of thymocytes. The results of the corresponding simulations are shown in Figs 2 and A.2. 4.1.2 Treatment modelling In this section, we will model the lymphoma treated with chemotherapy. The medications used are a combination of drugs, mainly Methotrexate and Purinethol. The treatment consists in several phases: induction, consolidation, re-induction and maintenance (Uyttebroeck et al., 2008). We model two treatment phases: the acute phase (induction) during 1–1.5 months and the maintenance phase, up to 2 years after the beginning of induction, (See Uyttebroeck et al. (2008)). The acute phase is not itself the focus of this work but its modelling, even simplified, is necessary. In the model, we consider a single cause of relapse: the presence of residual cancer cells that have not been killed during the acute treatment. It is therefore important to have patients with or without residual cancer cells when the maintenance treatment begins. The number of patient with potential relapses is a computation result and not a specified parameter in the model. Given the multiplicity of drugs, with different administration and doses, all the protocols of treatment cannot be realistically reproduced in the model. We consider a single drug summarizing the effect of all the drugs used. To allow this simplification, the periodicity of drug administration is also assumed to be sufficient to allow its concentration to be constant during the treatment. Cancer cells are exposed to the same quantity of drug regardless of their location. The maintenance treatment begins after the complete or partial disappearance of the tumour. We also assume that an equation similar to the equation used in the physiologically based pharmacokinetic model (PBPK) for Methotrexate and Mercaptopurine Ogungbenro (2014) can be used at the cellular level:   dp∗dt =k+∗p0∗−k−∗p∗. (4.1) here $$p^*$$ is the intracellular drug concentration and $$p_0^*$$ is the drug concentration in the thymus, which will be considered constant. The coefficient $$k^{+*}$$ determines the ability of a cell to absorb the drug and $$k^{-*}$$ to degrade it. When $$p^*(t)$$ reaches the critical value $$p_{crit}^*$$, the cell dies (Fig. A.3). Cells have the ability to change their response to treatment by changing the value of $$p_{crit}^m$$. It corresponds to the appearance of drug resistance. All simulations are carried out with a constant value of $$p_{crit}^m$$, except for the simulations to study treatment resistance (see Section 4.1.4). The treatment dosages for Mercaptopurine (purinethol) is 50 $$mg/m^2$$ daily and for Methotrexate 20 $$mg/m^2$$ weekly. The intracellular drug concentration which is necessary to kill malignant cells is not known. In what follows, both side of the equation (4.2) are divided by $$p_{crit}^m$$ in order to obtain dimensionless concentrations $$\frac{p_0^*}{p_{crit}^m}$$ and $$\frac{p^*}{p_{crit}^m}$$ , respectively notated $$p_0$$ and $$p$$. The value of $$p_{crit}^m$$ is set and the other parameters of the modelled chemotherapy are expressed as function of this value. The equation (4.2) becomes:   dpdt =k+p0−k−p. (4.2) Fig. A.4 shows the intracellular drug concentration versus time for various values of the parameters $$k^{-}$$ and $$k^{+}$$. The drug concentration is uniformly distributed inside the thymus and does not depend on time. Therefore, for a given simulation, all cells follow the same dynamics of drug accumulation. If the intracellular drug concentration reaches the critical level $$p_{crit}$$ at the end of the cell cycle, then the cell dies. Since cells differ by the duration of their cell cycle and by the value of their $$p_{crit}$$, some cells die while others survive even under treatment. $$p_{crit}$$ and $$T$$ are slightly different for all simulated cancer cells. These differences at the cellular level result in the simulation of different patient behaviours. The effect of chemotherapy on the thymus is illustrated in Fig. 3. Malignant cells are gradually eliminated. New precursors enter the thymus and recreate its normal structure. In Fig. 3, the drug dose $$p_0$$ was sufficient to kill all the cancer cells after a single cellular division. During computation, the screenshot increases dramatically the duration of simulations. This explains why we chose a drug dose, much higher than the usually administered dose, in order to show a quick disappearance of cancer cells and appearance of newly healthy cells. However, this simulation is not taken into account for the evaluation of the number of relapses (shown in Fig. 6) in order to get realistic results in accordance with the drug dosage in practice. 4.1.3 Simulation of patient’s response to treatment In the model after the acute phase therapy, the tumour regresses partially (residual disease) or totally (complete healing). When the maintenance treatment is administered, the number of cancer cells can be small or equal to zero. All the residual cancer cells may be killed by the maintenance treatment which will lead to a complete healing or they may survive to the maintenance treatment which will lead to a relapse. In our model, a tumour caused by cancer cells other than residual cancer cells is a new disease and not a relapse. This assumption seems to be close to the medical reality. Figure 4 illustrates the ability of the virtual patients to be cured or to relapse. The number of cancer cells is observed for a period almost equal to 2 years. The treatment regimens modelled are the same for two pairs of a virtual couple of patients, 24 months for the couple of patients A, 12 months for the couple of patients B. During the acute phase of the disease, the number of cancer cells increases. For patients A and B and for the two pairs, 1 month after the beginning of the treatment the number of cancer cells decreases and approaches zero. Figures 4 (a) and (c) show a pair of cured patients, Figs 4 (b) and (d) show one patient completely cured and one patient who relapsed. The relapse of patient B of the pair 1 is due to residual malignant cells that have not been killed during the maintenance treatment. This is due to the fact that cancer cells with a high value of $$p_{crit}$$ are predominant in this patient. Fig. 4. View largeDownload slide Numerical simulations of lymphoma therapy for two pairs of patients. The first vertical line shows the beginning of treatment of the acute phase, the second vertical line shows the end of the acute treatment and the beginning of the maintenance treatment, the third vertical line shows the end of the maintenance treatment for the patients in Figs (b) and (d). In those examples of patient simulations, we decided to stop treatment after 12 months for (b) and (d) patients. Other durations have been studied and results are shown in Fig. 6. Treatment is continued for patients in figs (a) and (c). Simulations made with the basic parameters values. Fig. 4. View largeDownload slide Numerical simulations of lymphoma therapy for two pairs of patients. The first vertical line shows the beginning of treatment of the acute phase, the second vertical line shows the end of the acute treatment and the beginning of the maintenance treatment, the third vertical line shows the end of the maintenance treatment for the patients in Figs (b) and (d). In those examples of patient simulations, we decided to stop treatment after 12 months for (b) and (d) patients. Other durations have been studied and results are shown in Fig. 6. Treatment is continued for patients in figs (a) and (c). Simulations made with the basic parameters values. 4.1.4 Resistance to treatment Cancer cells can develop a resistance to the treatment due to the selection and Darwinian evolution (Clairambault et al., 2013; Thomas et al., 2013). The most resistant cells will survive and proliferate while less resistant cells will die. In order to model this phenomenon, we introduce in the model the variation of the critical value $$p_{crit}^m$$. The value of $$p_{crit}^m$$ can be partly inherited from their mother cell, so that the value of $$p_{crit}^m$$ for daughter cell is equal to $$p_{crit} \pm \epsilon$$, with $$p_{crit}$$ the value of the mother cell. We set $$\epsilon=10\% p_{crit}^m$$. We have a variation of the critical value which leads to the gradual appearance of more resistant cells. Since less resistant cells will be eliminated by the treatment, the tumour will become more and more resistant to treatment. It is important to note that the initial drug concentration in newly born cells is set at zero because the drug half-life is very short. A specific simulation of 10 tumours, each arising from a single cancer cell allows to test the development of the cancer cells to resistance to the treatment (Fig. 5). Figure 5 shows the evolution of cell distribution with respect to their critical values of $$p_{crit}$$. We can conclude that these distributions move towards large critical values resulting in the appearance of cells resistant to treatment. Fig. 5. View largeDownload slide Cell distributions with respect to the values of $$p_{crit}$$ after 1 (blue), 2 (pink), 4 (yellow), 9 (brown), 14 (red), 20 (black) months. After several months of treatment, the mean value of $$p_{crit}$$ ($$p_{crit}^m$$) increases. We can observe a spreading of the Gaussian distribution, the standard deviation increases with time. Each histogram corresponds to an average of 10 simulations made with basic parameters values with the exception of $$p_{crit,daughter\text{ } cell}^m=p_{crit,mother\text{ }cell} \pm 10\% p_{crit,mother\text{ }cell}^m$$. (Color code refers to the electronic version of the paper.) Fig. 5. View largeDownload slide Cell distributions with respect to the values of $$p_{crit}$$ after 1 (blue), 2 (pink), 4 (yellow), 9 (brown), 14 (red), 20 (black) months. After several months of treatment, the mean value of $$p_{crit}$$ ($$p_{crit}^m$$) increases. We can observe a spreading of the Gaussian distribution, the standard deviation increases with time. Each histogram corresponds to an average of 10 simulations made with basic parameters values with the exception of $$p_{crit,daughter\text{ } cell}^m=p_{crit,mother\text{ }cell} \pm 10\% p_{crit,mother\text{ }cell}^m$$. (Color code refers to the electronic version of the paper.) 5. Simulation of virtual patients populations and comparison with patients from EuroLB database 5.1 Sensitivity of cancer cells to treatment Cancer relapse occurs due to remaining malignant cells after treatment. The efficacy of an anticancer drug depends on the administered dose and on the characteristics of the cell line. Quantitative parameters, such as growth inhibition ($$GI_{50}$$, the concentration for 50% of maximal inhibition of cell proliferation) for each cell line can be determined. In order to predict the response of cancer patients to chemotherapeutic agents, tumour chemosensitivity assay, drug response assays, or drug sensitivity assays are used (Ulukaya, 2006). Survival of cancer cells depends on the intracellular drug concentration (Durand and Olive, 1981). Cancer cell dies if it reaches some critical level $$p_{crit}$$. Biologically, the sensitivity depends on gene expression and, possibly, on some other factors. To our knowledge, the exact causes of relapses in LBL remains still unknown even if a genetic cause may be suspected. Patient response to treatment depends on the sensitivity of cancer cells. This assumption is consistent with the studies of in vivo and in vitro sensitivity to anticancer drugs. Pai et al. (2009) showed a positive correlation between the effect of methotrexate (one of the main drugs used for lymphoma treatment) on tumour cells in vitro and clinical response of the majority of patients. Similarly, the study of Konecny et al. (2000) reported a strong correlation between drug sensitivity/resistance and clinical response. They studied the in vitro chemosensitivity of primary epithelial ovarian cancer to drug combinations and compared the results with clinical response in patients. Similar tests can be used to determine in vitro chemosensitivity of lymphoma cells. Such tests would allow the determination of the critical value of the intracellular drug concentration $$p_{crit}$$ which is necessary to kill cancer cells. Since the data on the sensitivity of T-LBL cells to Methotrexate and Purinethol are not available, we consider the value of $$p_{crit}$$ as a parameter of the problem. We vary it and we model the response of the virtual population of patients to treatment assuming that patient response is determined by the sensitivity of cancer cells. 5.2 EuroLB database population The average duration of treatment of the 114 patients in the EuroLB database is 18 months. Short treatment durations (less than the mean value) are due to the death of patients, to modifications, to discontinuation of the treatment (e.g. because of side effects) or for unknown reasons. Among these real patients, 6 died from the disease. Data are available for three of them. They died at 216, 477 and 594 days (7, 16 and 20 months) after the beginning of the maintenance treatment. The corresponding data are summarized in Table 3. 5.3 Virtual populations Two populations of virtual patients have been constructed for this study. A first population of 122 pairs of patients allows to compare two durations of the maintenance treatment: 12 month versus 24 months. Table 1 shows the percentage of relapse of patients in the virtual population and of the patients in the EuroLB database population. Another population (two sets of 100 patients have been generated) allows to test the number of relapses as a function of the treatment durations and doses (data are summarized in Fig. 6). In the virtual population, we can see a relationship between the number of relapses and the duration of the treatment. In the simulations, for a duration of treatment less than 6 months, all the patients relapse. For a higher duration of treatment, the number of relapses decreases as a function of the duration of treatment and stabilizes around a low value (close to the value observed for a 2 years treatment). The duration necessary to obtain the stabilization of the number of relapses also depends on the choice of the parameters. There is a stabilization after 12 months of maintenance therapy for $$p_{0}=1.34$$ and after 18 months for $$p_{0}=1.3$$ (Fig. 6). Relationships and interactions between treatment duration, relapse rate and drug doses could be investigated in future studies. Table 1. Results of the simulation of a virtual population of 122 pairs of virtual patients (244 patients) and results observed in real patients (Uyttebroeck et al., 2008) Treatment stopped after:  12 months  24 months  Number of relapses  16  13  Number of recovery  106  109  Total  122  122  Rate of relapse (simulations)  13%  10%  Rate of relapse (in Uyttebroeck et al., (2008))  -  15%  Treatment stopped after:  12 months  24 months  Number of relapses  16  13  Number of recovery  106  109  Total  122  122  Rate of relapse (simulations)  13%  10%  Rate of relapse (in Uyttebroeck et al., (2008))  -  15%  Fig. 6. View largeDownload slide The proportion of relapses decreases (left) and the proportion of recoveries increases (right) as a function of the duration of the maintenance therapy for a constant $$p_{0}$$. Green curves correspond to $$p_{0}$$=1.34 (with a treatment of the acute phase equal to 1.25). Black curves correspond to $$p_{0}$$=1.3 (with a treatment of the acute phase equal to 1.3). Decrease of $$p_0$$ shifts the curves as follows: longer treatment is required to get the same ratio of relapsed and of recovered patients. Each point of the curve is the results of a group of 10 virtual patients, corresponding to a total of 100 patients per trial (left and right hand graphics: two populations of a hundred patients have been generated). Simulations made with basic parameters values with the exception of $$p_{0}$$. The proportion of relapses is not monotonically decreasing with the duration of the treatment because of random fluctuations determined by individual characteristics of the virtual patients. (Color code refers to the electronic version of the paper.) Fig. 6. View largeDownload slide The proportion of relapses decreases (left) and the proportion of recoveries increases (right) as a function of the duration of the maintenance therapy for a constant $$p_{0}$$. Green curves correspond to $$p_{0}$$=1.34 (with a treatment of the acute phase equal to 1.25). Black curves correspond to $$p_{0}$$=1.3 (with a treatment of the acute phase equal to 1.3). Decrease of $$p_0$$ shifts the curves as follows: longer treatment is required to get the same ratio of relapsed and of recovered patients. Each point of the curve is the results of a group of 10 virtual patients, corresponding to a total of 100 patients per trial (left and right hand graphics: two populations of a hundred patients have been generated). Simulations made with basic parameters values with the exception of $$p_{0}$$. The proportion of relapses is not monotonically decreasing with the duration of the treatment because of random fluctuations determined by individual characteristics of the virtual patients. (Color code refers to the electronic version of the paper.) 6. Discussion and perspectives Rare diseases are defined as a disease with a low prevalence that requires special combined efforts to address them. A global approach to rare diseases is consequently needed for establishing specific policies. There is a need to facilitate the development and availability of high quality, ethically researched, and appropriately authorized medicines in paediatric rare diseases, without subjecting children to unnecessary trials. The development of an in silico modelling approach and of a clinical trial simulation tool can help to choose the future most appropriate clinical trial(s) to evaluate therapeutic strategies. In the case of maintenance therapy for T-cell LBL, our mathematical model seems able to simulate the diversity of patient’s behaviour, complete healing or relapses that may occur during the treatment. In this work, we developed a model for T-LBL that could be further considered for phase III randomized clinical trials simulations purposes as considered in the CRESim project (Cornu et al., 2013; Nony et al., 2014) in order to compare trial design performances for rare diseases and rare events (relapses). Our mathematical model could be extended, considering for instance a time-varying quantity of drugs during a sequential therapy. After drug administration and depending on the pharmacokinetic properties of each drug, the intracellular drug concentration could be modelled and simulated. According to drug interactions and specific pharmacokinetic–pharmacodynamic relations, additional implications on the fate of cancer cells may occur. The simulation of sequential administrations of treatments, would therefore be closer to the real modalities of therapeutic protocols used in cancer patients. This could also be applied for the modelling and simulations of drug side effects. Another way to develop such a model could be further investigated, relapses being link to remaining cancer cells. The identification of biomarkers for the residual disease and their inclusion in the model would probably allow a more specific evaluation of the number of relapses. Generally, all simplifications of the disease and the therapeutic models may affect the results; accordingly the raw results of simulations cannot be directly applied for medical decisions in practice, for example the choice of duration of maintenance treatment. The purpose of our study was to build a model allowing to obtain results sufficiently close to the results observed in the study of Uyttebroeck et al. (2008) in terms of number of relapses. With the database, we know the value of the main output of the model (the % of relapses after a two years treatment) for real patients and we calibrate our model in order to achieve similar results. These results allow us to reduce the duration of virtual treatment. Let us also note that short treatments are not recommended because of the risks involved for patients. Since the probability of survival after a relapse in lymphoma is extremely low, then treatments with short durations are not used. According to our results, the simulation of therapeutic such as less than 6 months is clearly not recommended. Funding CRESim was funded by the Priomedchild Joint Call in 2010, REBUTTAL CRESim project no. 40-41800-98-008. The french part of the project was funded by ANR (agence nationale de la recherche). Members of the CRESim Project Group Leon Aarons, Agathe Bajard, Yves Bertrand, Frank Bretz, Daan Caudri, Charlotte Castellan, Sylvie Chabaud, Catherine Cornu, Frank Dufour, Nathalie Eymard, Roland Fisch, Renzo Guerrini, Behrouz Kassaï, Polina Kurbatova, Patrice Nony, Kayode Ogungbenro, David Pérol, Gérard Pons, Harm Tiddens, Anna Rosati. References Alarcón T. Byrneb H.M. & Mainia P.K. (2003) A cellular automaton model for tumour growth in inhomogeneous environment. J. Theor. Biol.,  225, 257– 274. Google Scholar CrossRef Search ADS   Anderson A.R.A. & Chaplain M.A.J. (1998) Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol.,  60, 857– 99. Google Scholar CrossRef Search ADS   Anderson A.R.A., Chaplain M.A.J. & Rejniak K.A. (2007) The cellular Potts model and its variants. Single-Cell-Based Models in Biology and Medicine. Basel, Boston, Berlin: Birkhäuser. ISBN 978-3-7643-8123-813. Bessonov N. Demin I. Pujo-Menjouet L. & Volpert V. (2009) A multi-agent model describing self-renewal of differentiation effects on the blood cell population. Math. Comput. Model.,  49, 2116– 2127. Google Scholar CrossRef Search ADS   Bessonov N. Crauste F. Fischer S. Kurbatova P. & Volpert V. (2011) Application of Hybrid Models to Blood Cell Production in the Bone Marrow. Math. Model. Nat. Phenom,  6, 2– 12. Google Scholar CrossRef Search ADS   Boehm T. (2012) Self-renewal of thymocytes in the absence of competitive precursor replenishment. J. Exp. Med.,  209, 1397– 1400. Google Scholar CrossRef Search ADS   Brousse N. (2009) Lymphome Lymphoblastique T enfant et adulte . Oral Communication. Brugières L. Minard V. & Patte C. (2012) Lymphomas in children and adolescents. Rev. Prat.,  62, 453– 8. Chiaretti S. & Foà R. (2009) T-cell acute lymphoblastic leukemia. Haematologica,  94, 160– 162. Google Scholar CrossRef Search ADS   Chrobak J.M. Bodnar M. & Herrero H. (2012) About a generalized model of lymphoma. J. Math. Anal. Appl.,  386, 813– 829. Google Scholar CrossRef Search ADS   Clairambault J. Magal P. & Volpert V. (2013) Cancer as evolutionary process. Eur Commun. Math. Theor. Biol.,  16, 19– 23. Clapp G. & Levy D. (2015) A review of mathematical models for leukemia and lymphoma. 116, 1– 6. Cornu C. Kassai B. Fisch R. Chiron C. Alberti C. Guerrini R. Rosati A. Pons G. Tiddens H. Chabaud S. Caudri D. Ballot C. Kurbatova P. Castellan A. Bajard A. Nony P.; CRESim & Epi-CRESim Project Groups. (2013) Experimental designs for small randomised clinical trials: an algorithm for choice. Orphanet. J. Rare Dis.,  8, 48. Google Scholar CrossRef Search ADS   Cortelazzo S. Ponzoni M. Ferreri A.J. & Hoelzer D. (2011) Lymphoblastic lymphoma. Crit. Rev. Oncol. Hematol.,  79, 330– 343. Google Scholar CrossRef Search ADS   Crompton T. Gilmour K.C. & Owen M.J. (1996) The MAP Kinase Pathway Controls Differentiation from Double-Negative to Double-Positive Thymocyte. Cell,  86, 243– 251. Google Scholar CrossRef Search ADS   Dallon J.C. (2007) Models with lattice-free center-based cells interacting with continuum environment variables. Single-cell-based models in biology and medecine , 197– 219. Deutsch A. (2007) Lattice-gas cellular automaton modelling of morphogenetic motion in developing cell systems In: Single Cell Based Models in Biology and Medicine, Basel: Birkhäuser, 29– 51. Durand R.E. & Olive P.L. (1981) Flow Cytometry Studies of Intracellular Adriamycin in Single Cells in single cells in vitro. Cancer Res. , 3489– 3494. Eymard N. Bessonov N. Gandrillon O. Koury M.J. & Volpert V. (2015) The role of spatial organization of cells in erythropoiesis. J. Math. Biol.,  70, 71– 97. Google Scholar CrossRef Search ADS   Eymard N. & Kurbatova P. (2015) Hybrid models in hematopoiesis. Math. Model. Nat. Phenom.,  10, 48– 63. Google Scholar CrossRef Search ADS   Fischer S. Kurbatova P. Bessonov N. Gandrillon O. Volpert V. & Crauste F. (2012) Modelling erythroblastic islands: using a hybrid model to assess the function of central macrophage. J. Theor. Biol.,  298, 92– 106. Google Scholar CrossRef Search ADS   Frieboes H.B., Smith B.R Chuang Yao-Li Ito K. Roettgers A.M Gambhir S,S. & Cristini V. (2013) An Integrated Computational/Experimental Model of Lymphoma Growth. PLOS Comput. Biol. , 9, e1003008. Hochberg J. & Cairo M.S. (2009) Childhood and adolescent lymphoblastic lymphoma: end of the beginning and future directions. Pediatr. Blood Cancer,  53, 917– 919. Google Scholar CrossRef Search ADS   Hoelzer D. & Gokbuget N. (2009) T-cell lymphoblastic lymphoma and T-cell acute lymphoblastic leukemia: a separate entity? Clin. Lymphoma Myeloma.,  9 (Suppl. (3)), S 214– 221. Google Scholar CrossRef Search ADS   Janeway C.A. et al.   (2011) Immunobiology, 6 th edition. Taylor & Francis Group, UK: Garland Science. Konecny G. Crohns C. Pegram M. Felber M. Lude S. Kurbacher C. Cree I.A Hepp H. & Untch M. (2000) Correlation of drug response with the ATP tumorchemosensitivity assay in primary FIGO stage III ovarian cancer. Gynecol. Oncol.,  77, 258– 263. Google Scholar CrossRef Search ADS   Kurbatova P. Bernard S. Bessonov N. Crauste F. Demin I. Dumontet C. Fischer S. & Volpert V. (2011) Hybrid model of erythropoiesis and leukemia treatment with cytosine arabinoside. SIAM J. App. Math.,  71 (6), 2246– 2268. Google Scholar CrossRef Search ADS   Kurbatova P. Eymard N. & Volpert V. (2013) Hybrid model of erythropoiesis. Acta Biotheoretica , 61, 305– 315. Google Scholar CrossRef Search ADS PubMed  Ma F. Manabe A. Wang D. Ito M. Kikuchi A. Wada M. Ito M. Ohara A. Hosoya R. Asano K. & Tsuji K. (2002) Growth of human T cell acute lymphoblastic leukemia lymphoblasts in NOD/SCID mouse fetal thymus organ culture. Leukemia,  16, 1541– 1548. Google Scholar CrossRef Search ADS   Martins V.C Ruggiero E. Schlenner S.M Madan V. Schmidt M. Fink P.J von Kalle C. & Rodewald H.-R. (2012) Thymus-autonomous T cell development in the absence of progenitor import. J. Exp. Med.,  209, 1409– 1417. Google Scholar CrossRef Search ADS   Nony P. Kurbatova K. Bajard A. Malik S. Castellan C. Chabaud S. Volpert V. Eymard N. Kassai B. Cornu C. & The CRESim and Epi-CRESim study groups. (2014) A methodological framework for drug development in rare diseases. Orphanet J. Rare Dis.,  9, 164. Google Scholar CrossRef Search ADS   Ogungbenro K. & Aarons L. (2014) Physiologically based pharmacokinetic modelling of Methotrexate and 6-Mercaptopurine in adults and children. J. Pharmacokinet. Pharmacodyn.,  41, 173– 185. Google Scholar CrossRef Search ADS   Page K.M & Uhr J.W. (2005) Mathematical models of cancer dormancy. Leuk. Lymphoma,  46, 313– 327. Google Scholar CrossRef Search ADS   Pai R.B Lalitha R.M Pai S.B Kumaraswamy S.V Lalitha N. & Bhargava M.K. (2009) Analysis of in vitro and vivo sensibility of oral cancer cells to methotrexate. Exp. Oncol.,  31, 118– 120. Peaudecerf L. Lemos S. Galgano A. Krenn G. Vasseur F. Di Santo J.P Ezine S. & Rocha B. (2012) Thymocytes may persist and differentiate without any input from bone marrow progenitors. J. Exp. Med.,  209, 1401– 1408. Google Scholar CrossRef Search ADS   Ribba B. Marron K. Agur Z. Alarc&#x2018;on T. & Maini P.K. (2005) A mathematical model of Doxorubicin treatment efficacy for non-Hodgkin’s lymphoma: investigation of the current protocol through theoretical modelling results. Bull. Math. Biol.,  67, 79– 99. Google Scholar CrossRef Search ADS   Rincon M. Flavell R.A. & Davis R.J. (2001) Signal transduction by MAP kinases in T lymphocytes. Oncogene,  20, 2490– 2497. Google Scholar CrossRef Search ADS   Schmidt E. & Burkhardt B. (2013) Lymphoblastic lymphoma in childhood and adolescence. Pediatr. Hematol. Oncol.,  30, 484– 508. Google Scholar CrossRef Search ADS   Thomas F. Fisher D. Fort P. Marie J-P Daoust S. Roche B. Gruna C. Cosseau C. Mitta G. Baghdiguian S. Rousset F. Lassus P. Assenat E. Grégoire D. Missé D. Lorz A. Billy F. Vainchenker W. Delhommeau F. Koscielny S. Itzykson, R., Tang R. Fava F. Ballesta A. Lepoutre T. Krasinska L. Dulic V. Raynaud P. Blache P. Quittau-Prevostel C. Vignal E. Trauchessec H. Perthame B. Clairambault J. Volpert V. Solary E. Hibner U. & Hochberg M.E. (2013) Applying ecological and evolutionary theory to cancer: a long and winding road. Evol. Appl.,  6, 1– 10. ISSN 1752– 4571. Google Scholar CrossRef Search ADS   Ulukaya E. (2006) Drug response assay: an increasing trend in designation of tailored-chemotherapy for more rational management of cancer patients. Adv. Mol. Med ., 2, 53– 58. Uyttebroeck A. Suciu S. Laureys G. Robert A. Pacquement H. Ferster A. Marguerite G. Mazingue F. Renard M. Lutz P. Rialland X. Mechinaud F. Cavé H. Baila L. & Bertrand Y. (2008) Treatment of childhood T-cell lymphoblastic lymphoma according to the strategy for acute lymphoblastic leukaemia, without radiotherapy: long term results of the EORTC CLG 58881 trial. Eur. J. Cancer,  44, 840– 846. Google Scholar CrossRef Search ADS   Appendix A A.1 Implementation of numerical algorithms The software used for the simulation of lymphoma is written in $$C^{++}$$, with operating system Ubuntu. The computations are carried out in dimensionless length units in such a way that the initial cell diameter corresponds to one unit. The computation domain is a square with the side equal 100 length units. In dimensional variables we consider cell diameters to be 10 microns. The structure of the software is as follow: Creation of a cell culture (list of cells) and of the computation domain. For each cell: Initialization of parameters. Birth of the healthy cells of a section of thymus and of one cancer cell. Each cell has a duration of life, coordinates, and intracellular concentrations. Forces between cells and coordinate of cells are calculated, movement is due to Newton’s second law. Intracellular and extracellular concentrations are calculated. Application of the rules governing healthy cell divisions. Every time step, age of cell is evaluated. While the age of the cell is less than or equal to the cell cycle, the cell grows. If a cell survives to the three selections, it self-renews or it differentiates according to the intracellular concentration. Application of the rules governing cancer cell divisions If the age of a cancer cell is equal or exceeds cell cycle duration and if the quantity of treatment is less than a threshold, the cell divides. Otherwise the cell dies. Treatment of the disease based of the date. Updated of the cell culture. A.2 Values of parameters Three kinds of parameters are used for this study: Individual-based modelling creates many implicit parameters, cell radius, cell mass, parameters of motion, etc. that are not essential in the results of simulations. They are related to the graphical outputs of the implemented model. Their values are expressed in arbitrary units. Some parameters clearly change the overall cell behaviour. In our model, a single protein $$u$$ and a single extracellular substance $$G$$ are considered as a way to obtain the three behaviours of a cell: apoptosis, self renewal and differentiation. Their concentrations are described by ODE and PDE but there are no available experimental data to set specific numerical values of parameters in these equations. Only expected division rates are accessible. In this case, we used dichotomy algorithms to obtain parameter values allowing division rates close to reality. On the other hand, the diffusion coefficient for growth factors in the reaction-diffusion equation is chosen in order to locate growth factors near the nurse cells. Moreover, the range of diffusion is set according to the size of the cells in order to obtain a realistic process. The modelled treatment is a combination of drugs with a mean value of toxicity for cancer cells (denoted $$p_{crit}^m$$ ) that cannot be exactly determined and that varies from one cell to another. For this reason, $$p_{crit}^m$$ is an arbitrary value and the intracellular concentration of treatment is normalized by this value. Therefore, the lack of knowledge of the lethal value for each cell does not affect the simulation process. For similar reasons (to our knowledge, it seems impossible to estimate the drung intake and metabolism in each cancer cell) the coefficients $$k^+$$ and $$k^{-}$$ are also function of the value $$p_{crit}^m$$. Table A.1 gives the basic value of parameters used for all the presented simulations unless otherwise specified in their legends. The unit of the dimensionless concentration is denoted $$C$$, dimensionless mass unit is denoted by M and L is dimensionless lengh unit. Table A.1 Basic values of parameters Parameter  Value  Unit  Cell cycle length  72  h (hours)  Cell cycle variation  9  h (hours)  $$r_0$$  0.01  $$L$$  $$m$$  1  $$M$$  $$\mu$$  $$6 \times 10^5$$  $$h^{-1}$$  initial value of $$u$$  0.01  $$C$$  $$c_{0}$$ (healthy cells)  $$6.905 \times 10^{-6}$$  $$h^{-1}$$  $$c_{0}$$ (cancer cells)  $$6.905 \times 10^{-4}$$  $$h^{-1}$$  $$k$$  0.005  $$C^{-1}$$  $$u_{crit}$$  1  $$C$$  $$k^{+*}$$  0.001  $$C^{-1}$$  $$k^{-*}$$  0.001  $$C^{-1}$$  $$p_{crit}^m$$  3  $$C$$  $$p_0^*$$ acute phase  3.75  $$C$$  $$p_0^*$$ maintenance phase  4  $$C$$  $$\gamma$$  0.005  $$h^{-1}$$  $$D$$  $$0.5 \times 10^{-4}$$  $$L^2.h^{-1}$$  $$w$$  $$0.5 \times 10^{-3}$$  molecules / $$L^2$$ / $$h$$  Parameter  Value  Unit  Cell cycle length  72  h (hours)  Cell cycle variation  9  h (hours)  $$r_0$$  0.01  $$L$$  $$m$$  1  $$M$$  $$\mu$$  $$6 \times 10^5$$  $$h^{-1}$$  initial value of $$u$$  0.01  $$C$$  $$c_{0}$$ (healthy cells)  $$6.905 \times 10^{-6}$$  $$h^{-1}$$  $$c_{0}$$ (cancer cells)  $$6.905 \times 10^{-4}$$  $$h^{-1}$$  $$k$$  0.005  $$C^{-1}$$  $$u_{crit}$$  1  $$C$$  $$k^{+*}$$  0.001  $$C^{-1}$$  $$k^{-*}$$  0.001  $$C^{-1}$$  $$p_{crit}^m$$  3  $$C$$  $$p_0^*$$ acute phase  3.75  $$C$$  $$p_0^*$$ maintenance phase  4  $$C$$  $$\gamma$$  0.005  $$h^{-1}$$  $$D$$  $$0.5 \times 10^{-4}$$  $$L^2.h^{-1}$$  $$w$$  $$0.5 \times 10^{-3}$$  molecules / $$L^2$$ / $$h$$  A.3 Results of treatment Table A.2 Results from EuroLB database (follow up of 144 patients during two years Duration of treatment  24 months  Number of deaths  6 ( 5% )  Number of patients  114  Duration of treatment  24 months  Number of deaths  6 ( 5% )  Number of patients  114  A.4 Additional figures Figures A.1, A.2 and A.4 show simulations using basic parameters values (Table A.1). Fig. A.1. View largeDownload slide Numerical simulations describing a healthy section of the thymus cortex. The number of cells (DN (blue), DP (green), SP (red)) in a section of healthy thymus during one year, in the computational domain. (Color code refers to the electronic version of the paper.) Fig. A.1. View largeDownload slide Numerical simulations describing a healthy section of the thymus cortex. The number of cells (DN (blue), DP (green), SP (red)) in a section of healthy thymus during one year, in the computational domain. (Color code refers to the electronic version of the paper.) Fig. A.2. View largeDownload slide The number of cells before and during the development of a tumour, while the number of malignant cells grows exponentially, the number of normal cells decreases. The curve showing the number of malignant cells is in purple (malignant cells M), the corresponding curves for normal cells are in blue (DN), green (DP), red (simple positive SP) and black (precursors P). (Color code refers to the electronic version of the paper.) Fig. A.2. View largeDownload slide The number of cells before and during the development of a tumour, while the number of malignant cells grows exponentially, the number of normal cells decreases. The curve showing the number of malignant cells is in purple (malignant cells M), the corresponding curves for normal cells are in blue (DN), green (DP), red (simple positive SP) and black (precursors P). (Color code refers to the electronic version of the paper.) Fig. A.3. View largeDownload slide If the drug amount in a given cell exceeds $$p_{crit}^*$$, at the end of its life, the cell dies otherwise it divides. Fig. A.3. View largeDownload slide If the drug amount in a given cell exceeds $$p_{crit}^*$$, at the end of its life, the cell dies otherwise it divides. Fig. A.4. View largeDownload slide Ratio of intracellular drug concentration $$p(t)$$ and the threshold $$p_{crit}$$ during cell cycle for $$k_+=2 \times 10^{-3}$$, $$k_-$$=$$10^{-3}$$ (blue), $$k_+$$=$$k_-$$=$$10^{-3}$$ (black), $$k_+$$=$$10^{-3}$$, $$k_-=2 \times 10^{-3}$$ (green), ($$p_0=1.25$$). (Color code refers to the electronic version of the paper.) Fig. A.4. View largeDownload slide Ratio of intracellular drug concentration $$p(t)$$ and the threshold $$p_{crit}$$ during cell cycle for $$k_+=2 \times 10^{-3}$$, $$k_-$$=$$10^{-3}$$ (blue), $$k_+$$=$$k_-$$=$$10^{-3}$$ (black), $$k_+$$=$$10^{-3}$$, $$k_-=2 \times 10^{-3}$$ (green), ($$p_0=1.25$$). (Color code refers to the electronic version of the paper.) © The authors 2016. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical Medicine And Biology: A Journal Of The Ima Oxford University Press

# Mathematical model of T-cell lymphoblastic lymphoma: disease, treatment, cure or relapse of a virtual cohort of patients

, Volume 35 (1) – Mar 1, 2018
23 pages

/lp/ou_press/mathematical-model-of-t-cell-lymphoblastic-lymphoma-disease-treatment-Fa6b00K49C
Publisher
Oxford University Press
ISSN
1477-8599
eISSN
1477-8602
D.O.I.
10.1093/imammb/dqw019
Publisher site
See Article on Publisher Site

### Abstract

Abstract T lymphoblastic lymphoma (T-LBL) is a rare type of lymphoma with a good prognosis with a remission rate of 85%. Patients can be completely cured or can relapse during or after a 2-year treatment. Relapses usually occur early after the remission of the acute phase. The median time of relapse is equal to 1 year, after the occurrence of complete remission (range 0.2–5.9 years) (Uyttebroeck et al., 2008). It can be assumed that patients may be treated longer than necessary with undue toxicity. The aim of our model was to investigate whether the duration of the maintenance therapy could be reduced without increasing the risk of relapses and to determine the minimum treatment duration that could be tested in a future clinical trial. We developed a mathematical model of virtual patients with T-LBL in order to obtain a proportion of virtual relapses close to the one observed in the real population of patients from the EuroLB database. Our simulations reproduced a 2-year follow-up required to study the onset of the disease, the treatment of the acute phase and the maintenance treatment phase. 1. Introduction This work was carried out as a part of the CRESim study (Child-Rare-Euro-Simulation), a research project funded by the ERA-NET PrioMedChild (Priority Medicines for Children). The main objective of the CRESim project was to create a platform for performing in silico experiments assessing and optimizing designs of RCTs for drug evaluation in children with rare diseases. Virtual RCTs can be simulated using different experimental designs and compared in terms of trial duration (for the patient, the investigator and the sponsor) and precision of the estimation of the treatment effect. In this particular work, the CRESim methodology was applied to study the duration of the maintenance treatment of T-cell lymphoblastic lymphoma. Lymphoblastic lymphoma (LBL) is a rare form of neoplasm that can develop from immature precursor cells in the bone marrow (BM) or from thymic T cells at varying stages of differentiation (Cortelazzo et al. (2011); Schmidt & Burkhardt (2013)). In developed countries, the annual cancer incidence in children under 15 years of age is estimated at 140 new cases per million and the LBL annual incidence at 0.3 to 0.5 cases per 100,000 children and adolescents (Brugiéres et al. (2012); Hochberg & Cairo (2009)). T lymphoblastic lymphoma (T-LBL) is a rare disease with few mathematical models developed for this type of cancer. A competition between normal and tumourous clonotypes of naive T-cells has been modelled in Chrobak et al. (2012) as a continuous-time bivariate Markov process. The study presented in Page & Uhr (2005) illustrates how the interaction between tumour and antibody for the murine BCL1 lymphoma leads to the dormancy of the tumour. The model takes into account cell cycle arrest and apoptosis of the tumour cells due to the immune response. The model of Non-Hodgkins lymphomas (NHL) (Ribba et al., 2005) includes a tumour structure incorporating mature and immature vessels, vascular structural adaptation and NHL cell-cycle kinetics in addition to pharmacokinetics and pharmacodynamics of the used drug (Doxorubicin). Another model of non-Hodgkin’ s lymphoma considers the tumour as a mixture of cells, interstitial fluid, and extracellular matrix (Frieboes et al., 2013). The concentrations of cells, as well as of oxygen and nutrients are described by partial differential equations. A recent review of lymphoma and leukaemia models (Clapp & Levy, 2014) focuses on haematopoiesis and drug improvement studies, and especially on drug resistance. Most T-LBL patients typically present with a mediastinal tumour. We therefore restrict our focus on a model of normal thymus, modelling three phases: its invasion by a tumour due to proliferation of a single cancer cell, the treatment of the acute phase and the maintenance treatment phase. For validation purposes, the results of our simulations compare the proportions of relapses in the virtual population with the ones from the EuroLB database of real patients. The aim is to build a model allowing the variation of the proportion of relapses of a population of T-LBL patients as a function of the drug dose and as a function of the treatment duration. The simulations should allow to test each of these essential elements separately or together. Inter-patients variability needed to reproduce different behaviours, cured or not, of real patients can be simulated by adding random noise to several existing models. However, to reach our objective, a more complex model is required. We use a mathematical model at a cellular level to simulate tumours and healthy thymuses; the presence or not, of surviving cancer cells after chemotherapy in a virtual patient determines his state of health. Given the clinical presentation of the disease, our model assumes a tumour arising in the thymus. We begin with the simulation of a healthy thymus. The model included the onset of the disease, the treatment of its acute and maintenance phase. In accordance with biomedical data, we assume that patient response to treatment depends on sensitivity of cancer cell to the intracellular drug concentration. In the first part of the work, we model the proliferation of a malignant tumour at a cell-scale. In this model, a given sensitivity of cancer cells to treatment is assigned to each cell as a parameter. Then, taking into account the probability density function of treatment sensitivity, the therapeutic effect of drug is simulated for virtual patients. 2. Thymus modelling The thymus is an important organ in the development of the immune system in young children. Its function becomes negligible in adults. The thymus is composed of several lobules, containing a cortex and a medulla separated by a cortico medullary junction. The thymus, in particular the cortex, is the site of maturation Brousse (2009), education and selection of T lymphocytes. Our simulations are focused on the processes called ontogenesis. 2.1 Healthy functioning thymus modelling: Ontogenesis The mathematical model of T cells ontogenesis takes into account four types of T cells: precursors, double negative (DN), double positive (DP) and simple positive (SP). The mechanisms allowing their progression between the consecutive steps implemented in our model are summarized in Fig. 1. The given percentages of apoptosis and the proportion of each kind of T cells are observed in reality and are computed by the model. In a healthy thymus, the proportion of T cells remains approximately constant: 3% of DN, 85% of DP and 12% of SP cells. Fig. 1. View largeDownload slide Discursive representation of the simulated ontogenesis. Fig. 1. View largeDownload slide Discursive representation of the simulated ontogenesis. Each day, $$5 \times 10^7$$ mature cells are produced but only $$10^6$$ cells leave the thymus. 98% of the thymocytes die by apoptosis by $$\beta$$, positive and negative selections. During their development, lymphocytes express several receptors allowing them to receive signals from their environment. Cells that fail to produce a functional pre-TCR (T Cell Membrane Receptor) are eliminated by apoptosis. This is called the $$\beta$$-selection. In positive selection, through contact with cortical epithelial cells, DN lymphocytes receive a signal of life or death. In negative selection, through contact with epithelial cells of medulla, DP surviving lymphocytes also receive this signal. The selection processes depend on the affinity TCRs for self MHC (Major Histocompatibility Complex) molecules from epithelial cells. Apoptosis occurs when there is a low affinity for positive selection and a high affinity for negative selection. Stem cells in the BM differentiate into thymic precursors and enter the thymus through blood vessels. After several stages of differentiation and proliferation, they then become T cells. Precursors are committed to a specific lineage through signalling. We only consider the T lineage by activation of their receptors $$Notch$$. During the early stages of maturation, precursors become DN1s, which are multipotent cells able to differentiate into T lymphocytes, NK (Natural Killer) lymphocytes and dendritic cells. Pro T thymocytes (DN2) mature and become DPs or NKs. Finally, DPs become SPs and exit the thymus through blood vessels (Janeway et al., 2011). This classification depends on the expression of their receptors CD4 and CD8, DP: $$CD4^+CD8^+$$, DN: $$CD4^-CD8^-$$, SP: $$CD4^+CD8^-$$ or $$CD4^-CD8^+$$. The maturation, migration, proliferation and differentiation of T cells lineage are controlled by interactions with other cells in the thymus, especially with cortical and medullary epithelial cells, nurse cells and dendritic cells. Our model only takes into account their actions on DNs, DPs and SPs and not their maturation processes. In the normal thymus, mature cells produce cytokines such as interleukins, growth factors and inhibitory factors. They provide a biological support for proliferation, differentiation or death signal. It is suggested that intracellular signals through different MAPK (mitogen-activated protein kinase) cascades selectively guide positive and negative selection of T lymphocytes (Hoelzer & Gokbuget, 2009). They are regulators of the differentiation of immature thymocytes from DN to DP cells, most probably acting as a transducer for pre-T cell receptor signalling (Crompton et al., 1996). The ERK (Extracellular Signal-regulated Kinease) and the p38 (a Mitogen-activated protein kinease) signalling cascades are also critically involved in TCR signals inducing positive selection and negative selection, respectively (Crompton et al., 1996). Differentiation, proliferation and self-renewal of real T cells are controlled by signals emerging from their microenvironment. In our model, the balance between cytokines families with antagonist roles is summarized by a reaction diffusion equation and by a death probability. At each differentiation stage a high rate of apoptosis is observed mostly due to positive and negative selections. Apoptosis of T cells, whatever the cause, is modelled by a probability of remaining alive for each cell type (cf. Subsection 3.1.4). Proliferation and differentiation of thymocytes are responses to signals sent by epithelial cells modelled by an extra regulation using partial differential equation (Subsections 3.1.2) that influences the T cells intraregulation simulated by an ordinary differential equation (cf. Subsection 3.1.3). 2.2 Disease modelling The disease begins in the thymus but the tumour and the effect of the treatment can also reach the whole body through blood vessels and affect the BM. The tumour results from a malignant transformation of clone abnormal progenitor cells with infinite self-renewal potential which occurs in the thymus (MA et al., 2002). Such cells present a genetic lesion that causes their proliferation and blocks their differentiation (Chiaretti & Foà R., 2009). Genetic aberrations in leukaemogenesis can affect transcription factor, protein-protein interaction, signal transduction, fate determination and cell cycle activator and differentiation (Chiaretti & Foà R., 2009). The consequences of these genetic mutations are very complex and are not modelled in our study. In 90% of cases, lymphomas are caused by an excessive proliferation of immature thymocytes Brousse (2009). Cancer cells appear to proliferate rather than to differentiate, causing the emergence and development of a tumour. Consequently, the tumour is modelled by a set of cells that proliferate with a low rate of apoptosis (probability of death close to zero) and without any differentiation. The intra and extraregulations, which induce differentiation, are not taken into account for the simulation of cancer cells. 3. Mathematical and computational model We visualize the simplified operating process of the thymus during two years using an hybrid model with an interface made with $$C^{++}$$ graphical library wxWidget (Fischer et al., 2012) included in specific software developed for haematopoiesis modelling (Bessonov et al., 2009). In this computed 2D representation, T cells are represented as circles with specific colours for each stage of differentiation, (Figs 2 and 3). Cell divisions are model as follow: during cell cycle, cell radius grows linearly. It becomes twice the initial radius at the end of the cycle. Then, each mother cell is replaced by two daughter cells. Geometrically, two circles with the initial radius are located inside the circle with a twice larger radius. The direction of cell division, that is the direction of the line connecting the centres of daughter cells is random (Fischer et al., 2012). Fig. 2. View largeDownload slide Development of the tumour in the thymus. The number of malignant cells grows exponentially with time. They invade the thymus and push out normal cells. The volume of the thymus increases and gradually cancer cells replace the normal T-cells. Cells shown in the figure are: DN (blue), DP (green), SP (red), stromal cells (large purple cells), malignant cells M (purple). The extracellular substances produced by nurse cells are shown as green halos. Simulations made with the basic parameters values. (Color code refers to the electronic version of the paper.) Fig. 2. View largeDownload slide Development of the tumour in the thymus. The number of malignant cells grows exponentially with time. They invade the thymus and push out normal cells. The volume of the thymus increases and gradually cancer cells replace the normal T-cells. Cells shown in the figure are: DN (blue), DP (green), SP (red), stromal cells (large purple cells), malignant cells M (purple). The extracellular substances produced by nurse cells are shown as green halos. Simulations made with the basic parameters values. (Color code refers to the electronic version of the paper.) Fig. 3. View largeDownload slide Simulations of lymphoma treatment by chemotherapy. Cancer cells are killed by the treatment. When the tumour disappears, precursors enter the thymus and recreate its normal structure after several stages of differentiation and proliferation. Precursors (gray), DN (blue), DP (green), SP (red), stromal cells (large purple cells) malignant cells (purple). Simulations made with basic parameters values with the exception of drug dose of acute treatment, equal to 5 (see explanation of this choice in the text). (Color code refers to the electronic version of the paper.) Fig. 3. View largeDownload slide Simulations of lymphoma treatment by chemotherapy. Cancer cells are killed by the treatment. When the tumour disappears, precursors enter the thymus and recreate its normal structure after several stages of differentiation and proliferation. Precursors (gray), DN (blue), DP (green), SP (red), stromal cells (large purple cells) malignant cells (purple). Simulations made with basic parameters values with the exception of drug dose of acute treatment, equal to 5 (see explanation of this choice in the text). (Color code refers to the electronic version of the paper.) On a computer view, a cell is a set of parameters (the meaning of all the considered parameters and their values are given in Appendix A.2). In order not to have a predetermine cell behaviour, these values can vary with time and cells location. Each cell has its own initial and critical values and coefficients of equations considered in the model. Their fate depend on the computation results done during their whole life (cf. Section 3.2). To increase legibility of this representation, cells cannot overlap. It is therefore necessary to consider their positions and displacements (Subsection 3.1.1, equations (3.4) and (3.3). A more complete approach of cells motion in hybrid models can be found in Kurbatova et al. (2013); Eymard et al. (2015); Eymard & Kurbatova (2015); Bessonov et al. (2011); Fischer et al. (2012). Displacements of cells in the thymus are important as it corresponds to a biologic reality; their maturation depends on their location in the thymus. The equations of the model are solved with an algorithm involving the spatial discretization of a computational domain, that is the considered thymus cut, and discretization of the time of simulation, that is each iteration of the computation corresponds to a time step $$\Delta t$$. In numerical simulations, equation (3.6) is solved by a finite difference method with Thomas algorithm and an alternative direction method. Neumann (no-flux) boundary conditions are considered as the boundary of the rectangular domain (Eymard et al., 2015; Fischer et al., 2012). Constant space and time steps are used. Equations (3.4) and (3.5) are solved by Euler method. The outputs of our model are the amounts of each kind of cells at any time. The results obtained with this approach cannot be directly compared with medical data. However, to our knowledge no medical examinations are able to provide such data. We defined the number of relapses at the end of the therapeutic maintenance phase simulation as follow: if the number of cancer cells is equal to zero, the patient is cured, otherwise the patient relapses. 3.1 T cell ontogenesis simulation 3.1.1 Hybrid discrete-continuous models Hybrid discrete-continuous models are widely used in the investigation of dynamics of cell populations in biological tissues and of organisms, which involve processes at different scales. This model is used to describe the education of T cells, the growth of the tumour and the action of the treatment on tumour cells. Hybrid models consist in coupling two models: one with a discrete approach for healthy and cancer cells modelled as objects interacting with their environment, and another with a continuous approach using ordinary and partial differential equations to describe intracellular and extracellular regulations that determine the fate of each cell (Alarcón et al., 2003; Anderson et al., 2007; Anderson & Chaplain, 1998; Dallon, 2007; Deutsch, 2007; Eymard et al., 2015; Kurbatova et al., 2013; Fischer et al., 2012). We consider the following equation to describe intracellular proteins in one cell:   dui(t)dt=F(ui(t),v(xi,t)), (3.1) where $$u_i$$ is a vector of intracellular concentrations of the cell $$i$$, $$v$$ is a vector of extracellular concentrations, $$F$$ is the vector of reaction rates which will be specified below and $$x_i$$ is the vector of coordinates of the centre of the cell. The quantity of protein depends on the concentrations of nutrients or biological products in the extracellular matrix described by the reaction-diffusion equation:   ∂v∂t=DΔv+H(v,c), (3.2) where $$c$$ is the local cell density, $$D$$ is the diffusion coefficient and $$H$$ is the rate of consumption or production of these substances by cells. These compounds can either be nutrients coming from outside or some other bio-chemical products consumed or produced by cells (Kurbatova et al., 2013; Eymard et al., 2015; Eymard & Kurbatova, 2015; Bessonov et al., 2011; Fischer et al., 2012; Kurbatova et al., 2011). The value of extracellular concentration $$v(x_i,t)$$ in equation (3.1) is taken at $$x=x_i$$. In our model, when cells divide, they push each other and cause their displacement inside the thymus. In order to describe mechanical interactions between cells, we restrict ourselves to the simplest model where cells are represented as compressible elastic spheres with an incompressible core (black circle inside each cell, Figs 2 and 3). Cell motion is determined by the sum of mechanical forces acting on a cell by other cells. Under the assumption of small deformations, the mechanical contact forces acting between cells can be expressed as a function $$f(d_{ij})$$ of the distance $$d_{ij}$$ between their centres. Cell motion is described by the displacement of its centre by Newton’s second law (Eymard & Kurbatova, 2015; Fischer et al., 2012):   mxi¨+μmxi˙−∑j≠if(dij)=0, (3.3) where $$m$$ is the mass of the cell, $$x_i$$ is the centre of the $$i$$-th cell, the second term in the left-hand side of this equation describes the friction by the surrounding medium, the last term represents the sum of potential forces acting on the i-th cell from all other cells (Bessonov et al., 2011; Fischer et al., 2012); (Kurbatova et al., 2011). The force between two spherical particles is given by the equation:   f(dij)={+∞,dij<d0−2H1κd0−dijdij−d0+2H1,dij<d0,0,dij≥d0  (3.4) where $$d_0$$ is the sum of cells radii, $$\kappa$$ is a positive parameter, and $$H_1$$ is the size of the outer compressible part of the cell. The force between the particles tends to infinity when $$d_{ij}$$ decreases to $$d_0-2H_1$$. The force is repulsive if $$d_{ij} < d_0$$ and is null otherwise. 3.1.2 Intracellular regulation of T cell In vivo experiments show that growth factors promote proliferation of early thymocyte cells (Martins et al., 2012; Peaudecerf et al., 2012). During their maturation, thymocyte cells move and away from nurse cells and they continue to differentiate. In order to describe the self-renewal and differentiation of thymocytes, we introduce an intracellular variable $$u$$ which corresponds to the concentration of the proteins MAPK. Its evolution is described by the equation:   dudt =c0+kG, (3.5) where $$c_0$$ is a positive constant and $$G$$ is the concentration of extracellular growth factors. If $$u$$ exceeds some critical value $$u_ {crit}$$ at the end of cell cycle, then the cell self-renews, otherwise it differentiates. 3.1.3 Extracellular regulation of T cell Early T cells are first located in a stromal cell niche and exposed to stromal factors (Boehm, 2012), which determine their proliferation and differentiation (Crompton et al., 1996; Rincon et al., 2001). Nurse cells produce growth factor $$G$$ which promotes self-renewal of thymocytes. If thymocyte cells are close to nurse cells, their intracellular concentration, denoted $$u$$, will increase due to the presence of growth factor. It will reach the critical value and the cell will self-renew. The cells located sufficiently far away from the nurse cells will differentiate. The concentration of the growth factor $$G$$ is described by the reaction-diffusion equation:   ∂G∂t=DΔG+W−γG, (3.6) where $$D$$ is the diffusion coefficient, $$W$$ is the production rate and $$\gamma$$ is the degradation rate blue of growth factor. These parameters are assumed to be constant. The values of parameters and initial conditions described in Subsections 3.1.2 and 3.1.3 are specified in Table 2 of Appendix A (Subsection A.2). This Table gives the basic value of parameters used for all the presented simulations unless otherwise specified in their legends. 3.1.4 $$\beta$$-selection, negative and positive selections Signal transductions are responsible for the different selections of T cells. Epithelial cells form a fibrous network in the thymus (Brousse, 2009). We assume that they are numerous enough and uniformly distributed in space such that each T cell can receive their signals of life and death from a cell to cell contact. In our simulations, apoptosis is modelled by a probability of remaining alive for each cell type. The probability of survival is equal to 60% for precursors, 40% for DN, 50% for DP and 40% for SP. After all maturational stages, the probability of 95% of apoptosis for healthy simulated cells is observed, this ratio is close to the ratio obtained by biological cells. 3.2 Thymus simulation The simulated cut of thymus is a set of cells, characterized by a list of parameters values. Initial values are assigned to cells at birth, such as: life expectancy $$T$$ (the term ‘cell cycle’ will also be used to define this duration), and initial concentration of protein or critical values such as lethal quantity of drug $$p_{crit}$$ and needed amount of intracellular protein for self renewal or differentiation $$u_{crit}$$. These parameters can be transmitted from mother to daughter cells or they can be random variables. At each iteration of the computation and for each cell, the position, the intracellular amount of drug $$p$$, of proteins $$u$$ and of growth factor available are updated, by solving equations (3.3, 3.4, 4.2, 3.5 and 3.6), respectively. It is important to note that most of the parameters are not biological characteristics of cells. They summarize the causes of a behaviour. Several proteins, growing factors and MAPK are involved in the choice of cell between apoptosis, self renewal and differentiation. In our model, a single protein $$u$$ and a single extracellular substance $$G$$ are considered. Similarly, the values of $$u_{crit}$$, $$k$$ and $$u_0$$ are chosen to obtain division rates close to the reality. The choice of the parameters’ values is a result of the calibration of the model by retrospective comparison of computational outputs with known real data. We compared the number of relapses after a two-year treatment observed in the real EuroLB patients and the virtual patients of our model. All other factors being equal, we varied the duration of treatments and the doses of therapy (cf. Section 4, results shown in Fig. 6). This can be considered as a sensitivity analysis of the most important inputs of the model: the duration of the treatment and its dose. However, a sensitivity analysis has not been performed for all others parameters, these simplifications are necessary to keep a reasonable computation time (the simulation of two patients requires about one day of the CPU time -bi-processeurs Intel Xeon 5650- of 10 nodes of computer cluster, namely more than 3 months of computations for all simulations). Cell cycle duration is taken in average of $$T=3$$ days plus/minus a random variation $$\delta$$. We set $$\delta = 1$$. Cell cycle varies in the time interval $$[T-\delta, T+\delta]$$, with a uniform distribution. After a period of time equal to its life expectancy, the cell chooses between self-renewal, differentiation or apoptosis. In the model, the amounts of intracellular proteins for healthy cells and of drug for cancer cells are compared with the corresponding threshold values and a decision is taken between the three options. When cells become SPs, they leave the thymus. To maintain a constant number of each cell types, new cells are introduced in the computational domain if there is some space available for them, by a periodic input of precursors. Consequently, new cells appear if the number of lymphocytes in the thymus is sufficiently low. In the case of tumour growth, the number of cancer cells is high and no new cells appear in the thymus. Spatial cell organization simulations with different cell types are shown in Fig. 2.Figure A.1 shows the results of numerical simulations of healthy thymus which correspond to one year of real time. The number of the different cell types is constant in average with some oscillations related to cell cycle. 4. Lymphoma development and treatment 4.1 Simulation of individual behaviour of a patient with a lymphoma: the onset of the disease and its treatment. 4.1.1 Simulation of spontaneous evolution of lymphoma The origin of a tumour is related to early thymocytes which have an excessive self-renewal potential and insufficient differentiation and apoptosis because of genetic mutations. To mimic tumour onset, we introduce a single mutated cell which proliferates. We model a tumour as a set of cancer cells. Malignant cells self-renew at the end of each cell cycle. They have a very low value of $$u_{crit}$$ and stay alive with a given probability $$q$$ ($$q=99\%$$ in all the simulations). In order to model the action of chemotherapy treatment on cancer cells, we prescribe to each cell the value $$p_{crit}$$ of the intracellular drug concentration needed to kill this particular cell. We characterize the cell population by an average value $$p_{crit}^m$$ with a normal distribution around this average value with standard deviation $$\sigma$$ equal 10% of the average value. In the computational domain, the multiplication of cancer cells increases the volume of the thymus and gradually invades it entirely. Healthy cells cannot develop normally and are progressively replaced by malignant cells. The tumour invades the thymus and destroys the spatial cell organization which is necessary for normal maturation of thymocytes. The results of the corresponding simulations are shown in Figs 2 and A.2. 4.1.2 Treatment modelling In this section, we will model the lymphoma treated with chemotherapy. The medications used are a combination of drugs, mainly Methotrexate and Purinethol. The treatment consists in several phases: induction, consolidation, re-induction and maintenance (Uyttebroeck et al., 2008). We model two treatment phases: the acute phase (induction) during 1–1.5 months and the maintenance phase, up to 2 years after the beginning of induction, (See Uyttebroeck et al. (2008)). The acute phase is not itself the focus of this work but its modelling, even simplified, is necessary. In the model, we consider a single cause of relapse: the presence of residual cancer cells that have not been killed during the acute treatment. It is therefore important to have patients with or without residual cancer cells when the maintenance treatment begins. The number of patient with potential relapses is a computation result and not a specified parameter in the model. Given the multiplicity of drugs, with different administration and doses, all the protocols of treatment cannot be realistically reproduced in the model. We consider a single drug summarizing the effect of all the drugs used. To allow this simplification, the periodicity of drug administration is also assumed to be sufficient to allow its concentration to be constant during the treatment. Cancer cells are exposed to the same quantity of drug regardless of their location. The maintenance treatment begins after the complete or partial disappearance of the tumour. We also assume that an equation similar to the equation used in the physiologically based pharmacokinetic model (PBPK) for Methotrexate and Mercaptopurine Ogungbenro (2014) can be used at the cellular level:   dp∗dt =k+∗p0∗−k−∗p∗. (4.1) here $$p^*$$ is the intracellular drug concentration and $$p_0^*$$ is the drug concentration in the thymus, which will be considered constant. The coefficient $$k^{+*}$$ determines the ability of a cell to absorb the drug and $$k^{-*}$$ to degrade it. When $$p^*(t)$$ reaches the critical value $$p_{crit}^*$$, the cell dies (Fig. A.3). Cells have the ability to change their response to treatment by changing the value of $$p_{crit}^m$$. It corresponds to the appearance of drug resistance. All simulations are carried out with a constant value of $$p_{crit}^m$$, except for the simulations to study treatment resistance (see Section 4.1.4). The treatment dosages for Mercaptopurine (purinethol) is 50 $$mg/m^2$$ daily and for Methotrexate 20 $$mg/m^2$$ weekly. The intracellular drug concentration which is necessary to kill malignant cells is not known. In what follows, both side of the equation (4.2) are divided by $$p_{crit}^m$$ in order to obtain dimensionless concentrations $$\frac{p_0^*}{p_{crit}^m}$$ and $$\frac{p^*}{p_{crit}^m}$$ , respectively notated $$p_0$$ and $$p$$. The value of $$p_{crit}^m$$ is set and the other parameters of the modelled chemotherapy are expressed as function of this value. The equation (4.2) becomes:   dpdt =k+p0−k−p. (4.2) Fig. A.4 shows the intracellular drug concentration versus time for various values of the parameters $$k^{-}$$ and $$k^{+}$$. The drug concentration is uniformly distributed inside the thymus and does not depend on time. Therefore, for a given simulation, all cells follow the same dynamics of drug accumulation. If the intracellular drug concentration reaches the critical level $$p_{crit}$$ at the end of the cell cycle, then the cell dies. Since cells differ by the duration of their cell cycle and by the value of their $$p_{crit}$$, some cells die while others survive even under treatment. $$p_{crit}$$ and $$T$$ are slightly different for all simulated cancer cells. These differences at the cellular level result in the simulation of different patient behaviours. The effect of chemotherapy on the thymus is illustrated in Fig. 3. Malignant cells are gradually eliminated. New precursors enter the thymus and recreate its normal structure. In Fig. 3, the drug dose $$p_0$$ was sufficient to kill all the cancer cells after a single cellular division. During computation, the screenshot increases dramatically the duration of simulations. This explains why we chose a drug dose, much higher than the usually administered dose, in order to show a quick disappearance of cancer cells and appearance of newly healthy cells. However, this simulation is not taken into account for the evaluation of the number of relapses (shown in Fig. 6) in order to get realistic results in accordance with the drug dosage in practice. 4.1.3 Simulation of patient’s response to treatment In the model after the acute phase therapy, the tumour regresses partially (residual disease) or totally (complete healing). When the maintenance treatment is administered, the number of cancer cells can be small or equal to zero. All the residual cancer cells may be killed by the maintenance treatment which will lead to a complete healing or they may survive to the maintenance treatment which will lead to a relapse. In our model, a tumour caused by cancer cells other than residual cancer cells is a new disease and not a relapse. This assumption seems to be close to the medical reality. Figure 4 illustrates the ability of the virtual patients to be cured or to relapse. The number of cancer cells is observed for a period almost equal to 2 years. The treatment regimens modelled are the same for two pairs of a virtual couple of patients, 24 months for the couple of patients A, 12 months for the couple of patients B. During the acute phase of the disease, the number of cancer cells increases. For patients A and B and for the two pairs, 1 month after the beginning of the treatment the number of cancer cells decreases and approaches zero. Figures 4 (a) and (c) show a pair of cured patients, Figs 4 (b) and (d) show one patient completely cured and one patient who relapsed. The relapse of patient B of the pair 1 is due to residual malignant cells that have not been killed during the maintenance treatment. This is due to the fact that cancer cells with a high value of $$p_{crit}$$ are predominant in this patient. Fig. 4. View largeDownload slide Numerical simulations of lymphoma therapy for two pairs of patients. The first vertical line shows the beginning of treatment of the acute phase, the second vertical line shows the end of the acute treatment and the beginning of the maintenance treatment, the third vertical line shows the end of the maintenance treatment for the patients in Figs (b) and (d). In those examples of patient simulations, we decided to stop treatment after 12 months for (b) and (d) patients. Other durations have been studied and results are shown in Fig. 6. Treatment is continued for patients in figs (a) and (c). Simulations made with the basic parameters values. Fig. 4. View largeDownload slide Numerical simulations of lymphoma therapy for two pairs of patients. The first vertical line shows the beginning of treatment of the acute phase, the second vertical line shows the end of the acute treatment and the beginning of the maintenance treatment, the third vertical line shows the end of the maintenance treatment for the patients in Figs (b) and (d). In those examples of patient simulations, we decided to stop treatment after 12 months for (b) and (d) patients. Other durations have been studied and results are shown in Fig. 6. Treatment is continued for patients in figs (a) and (c). Simulations made with the basic parameters values. 4.1.4 Resistance to treatment Cancer cells can develop a resistance to the treatment due to the selection and Darwinian evolution (Clairambault et al., 2013; Thomas et al., 2013). The most resistant cells will survive and proliferate while less resistant cells will die. In order to model this phenomenon, we introduce in the model the variation of the critical value $$p_{crit}^m$$. The value of $$p_{crit}^m$$ can be partly inherited from their mother cell, so that the value of $$p_{crit}^m$$ for daughter cell is equal to $$p_{crit} \pm \epsilon$$, with $$p_{crit}$$ the value of the mother cell. We set $$\epsilon=10\% p_{crit}^m$$. We have a variation of the critical value which leads to the gradual appearance of more resistant cells. Since less resistant cells will be eliminated by the treatment, the tumour will become more and more resistant to treatment. It is important to note that the initial drug concentration in newly born cells is set at zero because the drug half-life is very short. A specific simulation of 10 tumours, each arising from a single cancer cell allows to test the development of the cancer cells to resistance to the treatment (Fig. 5). Figure 5 shows the evolution of cell distribution with respect to their critical values of $$p_{crit}$$. We can conclude that these distributions move towards large critical values resulting in the appearance of cells resistant to treatment. Fig. 5. View largeDownload slide Cell distributions with respect to the values of $$p_{crit}$$ after 1 (blue), 2 (pink), 4 (yellow), 9 (brown), 14 (red), 20 (black) months. After several months of treatment, the mean value of $$p_{crit}$$ ($$p_{crit}^m$$) increases. We can observe a spreading of the Gaussian distribution, the standard deviation increases with time. Each histogram corresponds to an average of 10 simulations made with basic parameters values with the exception of $$p_{crit,daughter\text{ } cell}^m=p_{crit,mother\text{ }cell} \pm 10\% p_{crit,mother\text{ }cell}^m$$. (Color code refers to the electronic version of the paper.) Fig. 5. View largeDownload slide Cell distributions with respect to the values of $$p_{crit}$$ after 1 (blue), 2 (pink), 4 (yellow), 9 (brown), 14 (red), 20 (black) months. After several months of treatment, the mean value of $$p_{crit}$$ ($$p_{crit}^m$$) increases. We can observe a spreading of the Gaussian distribution, the standard deviation increases with time. Each histogram corresponds to an average of 10 simulations made with basic parameters values with the exception of $$p_{crit,daughter\text{ } cell}^m=p_{crit,mother\text{ }cell} \pm 10\% p_{crit,mother\text{ }cell}^m$$. (Color code refers to the electronic version of the paper.) 5. Simulation of virtual patients populations and comparison with patients from EuroLB database 5.1 Sensitivity of cancer cells to treatment Cancer relapse occurs due to remaining malignant cells after treatment. The efficacy of an anticancer drug depends on the administered dose and on the characteristics of the cell line. Quantitative parameters, such as growth inhibition ($$GI_{50}$$, the concentration for 50% of maximal inhibition of cell proliferation) for each cell line can be determined. In order to predict the response of cancer patients to chemotherapeutic agents, tumour chemosensitivity assay, drug response assays, or drug sensitivity assays are used (Ulukaya, 2006). Survival of cancer cells depends on the intracellular drug concentration (Durand and Olive, 1981). Cancer cell dies if it reaches some critical level $$p_{crit}$$. Biologically, the sensitivity depends on gene expression and, possibly, on some other factors. To our knowledge, the exact causes of relapses in LBL remains still unknown even if a genetic cause may be suspected. Patient response to treatment depends on the sensitivity of cancer cells. This assumption is consistent with the studies of in vivo and in vitro sensitivity to anticancer drugs. Pai et al. (2009) showed a positive correlation between the effect of methotrexate (one of the main drugs used for lymphoma treatment) on tumour cells in vitro and clinical response of the majority of patients. Similarly, the study of Konecny et al. (2000) reported a strong correlation between drug sensitivity/resistance and clinical response. They studied the in vitro chemosensitivity of primary epithelial ovarian cancer to drug combinations and compared the results with clinical response in patients. Similar tests can be used to determine in vitro chemosensitivity of lymphoma cells. Such tests would allow the determination of the critical value of the intracellular drug concentration $$p_{crit}$$ which is necessary to kill cancer cells. Since the data on the sensitivity of T-LBL cells to Methotrexate and Purinethol are not available, we consider the value of $$p_{crit}$$ as a parameter of the problem. We vary it and we model the response of the virtual population of patients to treatment assuming that patient response is determined by the sensitivity of cancer cells. 5.2 EuroLB database population The average duration of treatment of the 114 patients in the EuroLB database is 18 months. Short treatment durations (less than the mean value) are due to the death of patients, to modifications, to discontinuation of the treatment (e.g. because of side effects) or for unknown reasons. Among these real patients, 6 died from the disease. Data are available for three of them. They died at 216, 477 and 594 days (7, 16 and 20 months) after the beginning of the maintenance treatment. The corresponding data are summarized in Table 3. 5.3 Virtual populations Two populations of virtual patients have been constructed for this study. A first population of 122 pairs of patients allows to compare two durations of the maintenance treatment: 12 month versus 24 months. Table 1 shows the percentage of relapse of patients in the virtual population and of the patients in the EuroLB database population. Another population (two sets of 100 patients have been generated) allows to test the number of relapses as a function of the treatment durations and doses (data are summarized in Fig. 6). In the virtual population, we can see a relationship between the number of relapses and the duration of the treatment. In the simulations, for a duration of treatment less than 6 months, all the patients relapse. For a higher duration of treatment, the number of relapses decreases as a function of the duration of treatment and stabilizes around a low value (close to the value observed for a 2 years treatment). The duration necessary to obtain the stabilization of the number of relapses also depends on the choice of the parameters. There is a stabilization after 12 months of maintenance therapy for $$p_{0}=1.34$$ and after 18 months for $$p_{0}=1.3$$ (Fig. 6). Relationships and interactions between treatment duration, relapse rate and drug doses could be investigated in future studies. Table 1. Results of the simulation of a virtual population of 122 pairs of virtual patients (244 patients) and results observed in real patients (Uyttebroeck et al., 2008) Treatment stopped after:  12 months  24 months  Number of relapses  16  13  Number of recovery  106  109  Total  122  122  Rate of relapse (simulations)  13%  10%  Rate of relapse (in Uyttebroeck et al., (2008))  -  15%  Treatment stopped after:  12 months  24 months  Number of relapses  16  13  Number of recovery  106  109  Total  122  122  Rate of relapse (simulations)  13%  10%  Rate of relapse (in Uyttebroeck et al., (2008))  -  15%  Fig. 6. View largeDownload slide The proportion of relapses decreases (left) and the proportion of recoveries increases (right) as a function of the duration of the maintenance therapy for a constant $$p_{0}$$. Green curves correspond to $$p_{0}$$=1.34 (with a treatment of the acute phase equal to 1.25). Black curves correspond to $$p_{0}$$=1.3 (with a treatment of the acute phase equal to 1.3). Decrease of $$p_0$$ shifts the curves as follows: longer treatment is required to get the same ratio of relapsed and of recovered patients. Each point of the curve is the results of a group of 10 virtual patients, corresponding to a total of 100 patients per trial (left and right hand graphics: two populations of a hundred patients have been generated). Simulations made with basic parameters values with the exception of $$p_{0}$$. The proportion of relapses is not monotonically decreasing with the duration of the treatment because of random fluctuations determined by individual characteristics of the virtual patients. (Color code refers to the electronic version of the paper.) Fig. 6. View largeDownload slide The proportion of relapses decreases (left) and the proportion of recoveries increases (right) as a function of the duration of the maintenance therapy for a constant $$p_{0}$$. Green curves correspond to $$p_{0}$$=1.34 (with a treatment of the acute phase equal to 1.25). Black curves correspond to $$p_{0}$$=1.3 (with a treatment of the acute phase equal to 1.3). Decrease of $$p_0$$ shifts the curves as follows: longer treatment is required to get the same ratio of relapsed and of recovered patients. Each point of the curve is the results of a group of 10 virtual patients, corresponding to a total of 100 patients per trial (left and right hand graphics: two populations of a hundred patients have been generated). Simulations made with basic parameters values with the exception of $$p_{0}$$. The proportion of relapses is not monotonically decreasing with the duration of the treatment because of random fluctuations determined by individual characteristics of the virtual patients. (Color code refers to the electronic version of the paper.) 6. Discussion and perspectives Rare diseases are defined as a disease with a low prevalence that requires special combined efforts to address them. A global approach to rare diseases is consequently needed for establishing specific policies. There is a need to facilitate the development and availability of high quality, ethically researched, and appropriately authorized medicines in paediatric rare diseases, without subjecting children to unnecessary trials. The development of an in silico modelling approach and of a clinical trial simulation tool can help to choose the future most appropriate clinical trial(s) to evaluate therapeutic strategies. In the case of maintenance therapy for T-cell LBL, our mathematical model seems able to simulate the diversity of patient’s behaviour, complete healing or relapses that may occur during the treatment. In this work, we developed a model for T-LBL that could be further considered for phase III randomized clinical trials simulations purposes as considered in the CRESim project (Cornu et al., 2013; Nony et al., 2014) in order to compare trial design performances for rare diseases and rare events (relapses). Our mathematical model could be extended, considering for instance a time-varying quantity of drugs during a sequential therapy. After drug administration and depending on the pharmacokinetic properties of each drug, the intracellular drug concentration could be modelled and simulated. According to drug interactions and specific pharmacokinetic–pharmacodynamic relations, additional implications on the fate of cancer cells may occur. The simulation of sequential administrations of treatments, would therefore be closer to the real modalities of therapeutic protocols used in cancer patients. This could also be applied for the modelling and simulations of drug side effects. Another way to develop such a model could be further investigated, relapses being link to remaining cancer cells. The identification of biomarkers for the residual disease and their inclusion in the model would probably allow a more specific evaluation of the number of relapses. Generally, all simplifications of the disease and the therapeutic models may affect the results; accordingly the raw results of simulations cannot be directly applied for medical decisions in practice, for example the choice of duration of maintenance treatment. The purpose of our study was to build a model allowing to obtain results sufficiently close to the results observed in the study of Uyttebroeck et al. (2008) in terms of number of relapses. With the database, we know the value of the main output of the model (the % of relapses after a two years treatment) for real patients and we calibrate our model in order to achieve similar results. These results allow us to reduce the duration of virtual treatment. Let us also note that short treatments are not recommended because of the risks involved for patients. Since the probability of survival after a relapse in lymphoma is extremely low, then treatments with short durations are not used. According to our results, the simulation of therapeutic such as less than 6 months is clearly not recommended. Funding CRESim was funded by the Priomedchild Joint Call in 2010, REBUTTAL CRESim project no. 40-41800-98-008. The french part of the project was funded by ANR (agence nationale de la recherche). Members of the CRESim Project Group Leon Aarons, Agathe Bajard, Yves Bertrand, Frank Bretz, Daan Caudri, Charlotte Castellan, Sylvie Chabaud, Catherine Cornu, Frank Dufour, Nathalie Eymard, Roland Fisch, Renzo Guerrini, Behrouz Kassaï, Polina Kurbatova, Patrice Nony, Kayode Ogungbenro, David Pérol, Gérard Pons, Harm Tiddens, Anna Rosati. References Alarcón T. Byrneb H.M. & Mainia P.K. (2003) A cellular automaton model for tumour growth in inhomogeneous environment. J. Theor. Biol.,  225, 257– 274. Google Scholar CrossRef Search ADS   Anderson A.R.A. & Chaplain M.A.J. (1998) Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol.,  60, 857– 99. Google Scholar CrossRef Search ADS   Anderson A.R.A., Chaplain M.A.J. & Rejniak K.A. (2007) The cellular Potts model and its variants. Single-Cell-Based Models in Biology and Medicine. Basel, Boston, Berlin: Birkhäuser. ISBN 978-3-7643-8123-813. Bessonov N. Demin I. Pujo-Menjouet L. & Volpert V. (2009) A multi-agent model describing self-renewal of differentiation effects on the blood cell population. Math. Comput. Model.,  49, 2116– 2127. Google Scholar CrossRef Search ADS   Bessonov N. Crauste F. Fischer S. Kurbatova P. & Volpert V. (2011) Application of Hybrid Models to Blood Cell Production in the Bone Marrow. Math. Model. Nat. Phenom,  6, 2– 12. Google Scholar CrossRef Search ADS   Boehm T. (2012) Self-renewal of thymocytes in the absence of competitive precursor replenishment. J. Exp. Med.,  209, 1397– 1400. Google Scholar CrossRef Search ADS   Brousse N. (2009) Lymphome Lymphoblastique T enfant et adulte . Oral Communication. Brugières L. Minard V. & Patte C. (2012) Lymphomas in children and adolescents. Rev. Prat.,  62, 453– 8. Chiaretti S. & Foà R. (2009) T-cell acute lymphoblastic leukemia. Haematologica,  94, 160– 162. Google Scholar CrossRef Search ADS   Chrobak J.M. Bodnar M. & Herrero H. (2012) About a generalized model of lymphoma. J. Math. Anal. Appl.,  386, 813– 829. Google Scholar CrossRef Search ADS   Clairambault J. Magal P. & Volpert V. (2013) Cancer as evolutionary process. Eur Commun. Math. Theor. Biol.,  16, 19– 23. Clapp G. & Levy D. (2015) A review of mathematical models for leukemia and lymphoma. 116, 1– 6. Cornu C. Kassai B. Fisch R. Chiron C. Alberti C. Guerrini R. Rosati A. Pons G. Tiddens H. Chabaud S. Caudri D. Ballot C. Kurbatova P. Castellan A. Bajard A. Nony P.; CRESim & Epi-CRESim Project Groups. (2013) Experimental designs for small randomised clinical trials: an algorithm for choice. Orphanet. J. Rare Dis.,  8, 48. Google Scholar CrossRef Search ADS   Cortelazzo S. Ponzoni M. Ferreri A.J. & Hoelzer D. (2011) Lymphoblastic lymphoma. Crit. Rev. Oncol. Hematol.,  79, 330– 343. Google Scholar CrossRef Search ADS   Crompton T. Gilmour K.C. & Owen M.J. (1996) The MAP Kinase Pathway Controls Differentiation from Double-Negative to Double-Positive Thymocyte. Cell,  86, 243– 251. Google Scholar CrossRef Search ADS   Dallon J.C. (2007) Models with lattice-free center-based cells interacting with continuum environment variables. Single-cell-based models in biology and medecine , 197– 219. Deutsch A. (2007) Lattice-gas cellular automaton modelling of morphogenetic motion in developing cell systems In: Single Cell Based Models in Biology and Medicine, Basel: Birkhäuser, 29– 51. Durand R.E. & Olive P.L. (1981) Flow Cytometry Studies of Intracellular Adriamycin in Single Cells in single cells in vitro. Cancer Res. , 3489– 3494. Eymard N. Bessonov N. Gandrillon O. Koury M.J. & Volpert V. (2015) The role of spatial organization of cells in erythropoiesis. J. Math. Biol.,  70, 71– 97. Google Scholar CrossRef Search ADS   Eymard N. & Kurbatova P. (2015) Hybrid models in hematopoiesis. Math. Model. Nat. Phenom.,  10, 48– 63. Google Scholar CrossRef Search ADS   Fischer S. Kurbatova P. Bessonov N. Gandrillon O. Volpert V. & Crauste F. (2012) Modelling erythroblastic islands: using a hybrid model to assess the function of central macrophage. J. Theor. Biol.,  298, 92– 106. Google Scholar CrossRef Search ADS   Frieboes H.B., Smith B.R Chuang Yao-Li Ito K. Roettgers A.M Gambhir S,S. & Cristini V. (2013) An Integrated Computational/Experimental Model of Lymphoma Growth. PLOS Comput. Biol. , 9, e1003008. Hochberg J. & Cairo M.S. (2009) Childhood and adolescent lymphoblastic lymphoma: end of the beginning and future directions. Pediatr. Blood Cancer,  53, 917– 919. Google Scholar CrossRef Search ADS   Hoelzer D. & Gokbuget N. (2009) T-cell lymphoblastic lymphoma and T-cell acute lymphoblastic leukemia: a separate entity? Clin. Lymphoma Myeloma.,  9 (Suppl. (3)), S 214– 221. Google Scholar CrossRef Search ADS   Janeway C.A. et al.   (2011) Immunobiology, 6 th edition. Taylor & Francis Group, UK: Garland Science. Konecny G. Crohns C. Pegram M. Felber M. Lude S. Kurbacher C. Cree I.A Hepp H. & Untch M. (2000) Correlation of drug response with the ATP tumorchemosensitivity assay in primary FIGO stage III ovarian cancer. Gynecol. Oncol.,  77, 258– 263. Google Scholar CrossRef Search ADS   Kurbatova P. Bernard S. Bessonov N. Crauste F. Demin I. Dumontet C. Fischer S. & Volpert V. (2011) Hybrid model of erythropoiesis and leukemia treatment with cytosine arabinoside. SIAM J. App. Math.,  71 (6), 2246– 2268. Google Scholar CrossRef Search ADS   Kurbatova P. Eymard N. & Volpert V. (2013) Hybrid model of erythropoiesis. Acta Biotheoretica , 61, 305– 315. Google Scholar CrossRef Search ADS PubMed  Ma F. Manabe A. Wang D. Ito M. Kikuchi A. Wada M. Ito M. Ohara A. Hosoya R. Asano K. & Tsuji K. (2002) Growth of human T cell acute lymphoblastic leukemia lymphoblasts in NOD/SCID mouse fetal thymus organ culture. Leukemia,  16, 1541– 1548. Google Scholar CrossRef Search ADS   Martins V.C Ruggiero E. Schlenner S.M Madan V. Schmidt M. Fink P.J von Kalle C. & Rodewald H.-R. (2012) Thymus-autonomous T cell development in the absence of progenitor import. J. Exp. Med.,  209, 1409– 1417. Google Scholar CrossRef Search ADS   Nony P. Kurbatova K. Bajard A. Malik S. Castellan C. Chabaud S. Volpert V. Eymard N. Kassai B. Cornu C. & The CRESim and Epi-CRESim study groups. (2014) A methodological framework for drug development in rare diseases. Orphanet J. Rare Dis.,  9, 164. Google Scholar CrossRef Search ADS   Ogungbenro K. & Aarons L. (2014) Physiologically based pharmacokinetic modelling of Methotrexate and 6-Mercaptopurine in adults and children. J. Pharmacokinet. Pharmacodyn.,  41, 173– 185. Google Scholar CrossRef Search ADS   Page K.M & Uhr J.W. (2005) Mathematical models of cancer dormancy. Leuk. Lymphoma,  46, 313– 327. Google Scholar CrossRef Search ADS   Pai R.B Lalitha R.M Pai S.B Kumaraswamy S.V Lalitha N. & Bhargava M.K. (2009) Analysis of in vitro and vivo sensibility of oral cancer cells to methotrexate. Exp. Oncol.,  31, 118– 120. Peaudecerf L. Lemos S. Galgano A. Krenn G. Vasseur F. Di Santo J.P Ezine S. & Rocha B. (2012) Thymocytes may persist and differentiate without any input from bone marrow progenitors. J. Exp. Med.,  209, 1401– 1408. Google Scholar CrossRef Search ADS   Ribba B. Marron K. Agur Z. Alarc&#x2018;on T. & Maini P.K. (2005) A mathematical model of Doxorubicin treatment efficacy for non-Hodgkin’s lymphoma: investigation of the current protocol through theoretical modelling results. Bull. Math. Biol.,  67, 79– 99. Google Scholar CrossRef Search ADS   Rincon M. Flavell R.A. & Davis R.J. (2001) Signal transduction by MAP kinases in T lymphocytes. Oncogene,  20, 2490– 2497. Google Scholar CrossRef Search ADS   Schmidt E. & Burkhardt B. (2013) Lymphoblastic lymphoma in childhood and adolescence. Pediatr. Hematol. Oncol.,  30, 484– 508. Google Scholar CrossRef Search ADS   Thomas F. Fisher D. Fort P. Marie J-P Daoust S. Roche B. Gruna C. Cosseau C. Mitta G. Baghdiguian S. Rousset F. Lassus P. Assenat E. Grégoire D. Missé D. Lorz A. Billy F. Vainchenker W. Delhommeau F. Koscielny S. Itzykson, R., Tang R. Fava F. Ballesta A. Lepoutre T. Krasinska L. Dulic V. Raynaud P. Blache P. Quittau-Prevostel C. Vignal E. Trauchessec H. Perthame B. Clairambault J. Volpert V. Solary E. Hibner U. & Hochberg M.E. (2013) Applying ecological and evolutionary theory to cancer: a long and winding road. Evol. Appl.,  6, 1– 10. ISSN 1752– 4571. Google Scholar CrossRef Search ADS   Ulukaya E. (2006) Drug response assay: an increasing trend in designation of tailored-chemotherapy for more rational management of cancer patients. Adv. Mol. Med ., 2, 53– 58. Uyttebroeck A. Suciu S. Laureys G. Robert A. Pacquement H. Ferster A. Marguerite G. Mazingue F. Renard M. Lutz P. Rialland X. Mechinaud F. Cavé H. Baila L. & Bertrand Y. (2008) Treatment of childhood T-cell lymphoblastic lymphoma according to the strategy for acute lymphoblastic leukaemia, without radiotherapy: long term results of the EORTC CLG 58881 trial. Eur. J. Cancer,  44, 840– 846. Google Scholar CrossRef Search ADS   Appendix A A.1 Implementation of numerical algorithms The software used for the simulation of lymphoma is written in $$C^{++}$$, with operating system Ubuntu. The computations are carried out in dimensionless length units in such a way that the initial cell diameter corresponds to one unit. The computation domain is a square with the side equal 100 length units. In dimensional variables we consider cell diameters to be 10 microns. The structure of the software is as follow: Creation of a cell culture (list of cells) and of the computation domain. For each cell: Initialization of parameters. Birth of the healthy cells of a section of thymus and of one cancer cell. Each cell has a duration of life, coordinates, and intracellular concentrations. Forces between cells and coordinate of cells are calculated, movement is due to Newton’s second law. Intracellular and extracellular concentrations are calculated. Application of the rules governing healthy cell divisions. Every time step, age of cell is evaluated. While the age of the cell is less than or equal to the cell cycle, the cell grows. If a cell survives to the three selections, it self-renews or it differentiates according to the intracellular concentration. Application of the rules governing cancer cell divisions If the age of a cancer cell is equal or exceeds cell cycle duration and if the quantity of treatment is less than a threshold, the cell divides. Otherwise the cell dies. Treatment of the disease based of the date. Updated of the cell culture. A.2 Values of parameters Three kinds of parameters are used for this study: Individual-based modelling creates many implicit parameters, cell radius, cell mass, parameters of motion, etc. that are not essential in the results of simulations. They are related to the graphical outputs of the implemented model. Their values are expressed in arbitrary units. Some parameters clearly change the overall cell behaviour. In our model, a single protein $$u$$ and a single extracellular substance $$G$$ are considered as a way to obtain the three behaviours of a cell: apoptosis, self renewal and differentiation. Their concentrations are described by ODE and PDE but there are no available experimental data to set specific numerical values of parameters in these equations. Only expected division rates are accessible. In this case, we used dichotomy algorithms to obtain parameter values allowing division rates close to reality. On the other hand, the diffusion coefficient for growth factors in the reaction-diffusion equation is chosen in order to locate growth factors near the nurse cells. Moreover, the range of diffusion is set according to the size of the cells in order to obtain a realistic process. The modelled treatment is a combination of drugs with a mean value of toxicity for cancer cells (denoted $$p_{crit}^m$$ ) that cannot be exactly determined and that varies from one cell to another. For this reason, $$p_{crit}^m$$ is an arbitrary value and the intracellular concentration of treatment is normalized by this value. Therefore, the lack of knowledge of the lethal value for each cell does not affect the simulation process. For similar reasons (to our knowledge, it seems impossible to estimate the drung intake and metabolism in each cancer cell) the coefficients $$k^+$$ and $$k^{-}$$ are also function of the value $$p_{crit}^m$$. Table A.1 gives the basic value of parameters used for all the presented simulations unless otherwise specified in their legends. The unit of the dimensionless concentration is denoted $$C$$, dimensionless mass unit is denoted by M and L is dimensionless lengh unit. Table A.1 Basic values of parameters Parameter  Value  Unit  Cell cycle length  72  h (hours)  Cell cycle variation  9  h (hours)  $$r_0$$  0.01  $$L$$  $$m$$  1  $$M$$  $$\mu$$  $$6 \times 10^5$$  $$h^{-1}$$  initial value of $$u$$  0.01  $$C$$  $$c_{0}$$ (healthy cells)  $$6.905 \times 10^{-6}$$  $$h^{-1}$$  $$c_{0}$$ (cancer cells)  $$6.905 \times 10^{-4}$$  $$h^{-1}$$  $$k$$  0.005  $$C^{-1}$$  $$u_{crit}$$  1  $$C$$  $$k^{+*}$$  0.001  $$C^{-1}$$  $$k^{-*}$$  0.001  $$C^{-1}$$  $$p_{crit}^m$$  3  $$C$$  $$p_0^*$$ acute phase  3.75  $$C$$  $$p_0^*$$ maintenance phase  4  $$C$$  $$\gamma$$  0.005  $$h^{-1}$$  $$D$$  $$0.5 \times 10^{-4}$$  $$L^2.h^{-1}$$  $$w$$  $$0.5 \times 10^{-3}$$  molecules / $$L^2$$ / $$h$$  Parameter  Value  Unit  Cell cycle length  72  h (hours)  Cell cycle variation  9  h (hours)  $$r_0$$  0.01  $$L$$  $$m$$  1  $$M$$  $$\mu$$  $$6 \times 10^5$$  $$h^{-1}$$  initial value of $$u$$  0.01  $$C$$  $$c_{0}$$ (healthy cells)  $$6.905 \times 10^{-6}$$  $$h^{-1}$$  $$c_{0}$$ (cancer cells)  $$6.905 \times 10^{-4}$$  $$h^{-1}$$  $$k$$  0.005  $$C^{-1}$$  $$u_{crit}$$  1  $$C$$  $$k^{+*}$$  0.001  $$C^{-1}$$  $$k^{-*}$$  0.001  $$C^{-1}$$  $$p_{crit}^m$$  3  $$C$$  $$p_0^*$$ acute phase  3.75  $$C$$  $$p_0^*$$ maintenance phase  4  $$C$$  $$\gamma$$  0.005  $$h^{-1}$$  $$D$$  $$0.5 \times 10^{-4}$$  $$L^2.h^{-1}$$  $$w$$  $$0.5 \times 10^{-3}$$  molecules / $$L^2$$ / $$h$$  A.3 Results of treatment Table A.2 Results from EuroLB database (follow up of 144 patients during two years Duration of treatment  24 months  Number of deaths  6 ( 5% )  Number of patients  114  Duration of treatment  24 months  Number of deaths  6 ( 5% )  Number of patients  114  A.4 Additional figures Figures A.1, A.2 and A.4 show simulations using basic parameters values (Table A.1). Fig. A.1. View largeDownload slide Numerical simulations describing a healthy section of the thymus cortex. The number of cells (DN (blue), DP (green), SP (red)) in a section of healthy thymus during one year, in the computational domain. (Color code refers to the electronic version of the paper.) Fig. A.1. View largeDownload slide Numerical simulations describing a healthy section of the thymus cortex. The number of cells (DN (blue), DP (green), SP (red)) in a section of healthy thymus during one year, in the computational domain. (Color code refers to the electronic version of the paper.) Fig. A.2. View largeDownload slide The number of cells before and during the development of a tumour, while the number of malignant cells grows exponentially, the number of normal cells decreases. The curve showing the number of malignant cells is in purple (malignant cells M), the corresponding curves for normal cells are in blue (DN), green (DP), red (simple positive SP) and black (precursors P). (Color code refers to the electronic version of the paper.) Fig. A.2. View largeDownload slide The number of cells before and during the development of a tumour, while the number of malignant cells grows exponentially, the number of normal cells decreases. The curve showing the number of malignant cells is in purple (malignant cells M), the corresponding curves for normal cells are in blue (DN), green (DP), red (simple positive SP) and black (precursors P). (Color code refers to the electronic version of the paper.) Fig. A.3. View largeDownload slide If the drug amount in a given cell exceeds $$p_{crit}^*$$, at the end of its life, the cell dies otherwise it divides. Fig. A.3. View largeDownload slide If the drug amount in a given cell exceeds $$p_{crit}^*$$, at the end of its life, the cell dies otherwise it divides. Fig. A.4. View largeDownload slide Ratio of intracellular drug concentration $$p(t)$$ and the threshold $$p_{crit}$$ during cell cycle for $$k_+=2 \times 10^{-3}$$, $$k_-$$=$$10^{-3}$$ (blue), $$k_+$$=$$k_-$$=$$10^{-3}$$ (black), $$k_+$$=$$10^{-3}$$, $$k_-=2 \times 10^{-3}$$ (green), ($$p_0=1.25$$). (Color code refers to the electronic version of the paper.) Fig. A.4. View largeDownload slide Ratio of intracellular drug concentration $$p(t)$$ and the threshold $$p_{crit}$$ during cell cycle for $$k_+=2 \times 10^{-3}$$, $$k_-$$=$$10^{-3}$$ (blue), $$k_+$$=$$k_-$$=$$10^{-3}$$ (black), $$k_+$$=$$10^{-3}$$, $$k_-=2 \times 10^{-3}$$ (green), ($$p_0=1.25$$). (Color code refers to the electronic version of the paper.) © The authors 2016. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

### Journal

Mathematical Medicine And Biology: A Journal Of The ImaOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations