Mathematical modelling as a proof of concept for MPNs as a human inflammation model for cancer development

The chronic Philadelphia-negative myeloproliferative neoplasms (MPNs) are acquired stem cell neoplasms which ultimately may transform to acute myelogenous leukemia. Most recently, chronic inflammation has been described as an important factor for the development and progression of MPNs in the biological continuum from early cancer stage to the advanced myelofibrosis stage, the MPNs being described as “A Human Inflammation Model for Cancer Development“. This novel concept has been built upon clinical, experimental, genomic, immunological and not least epidemiological studies. Only a few studies have described the development of MPNs by mathematical models, and none have addressed the role of inflammation for clonal evolution and disease progression. Herein, we aim at using mathematical modelling to substantiate the concept of chronic inflammation as an important trigger and driver of MPNs.The basics of the model describe the proliferation from stem cells to mature cells including mutations of healthy stem cells to become malignant stem cells. We include a simple inflammatory coupling coping with cell death and affecting the basic model beneath. First, we describe the system without feedbacks or regulatory interactions. Next, we introduce inflammatory feedback into the system. Finally, we include other feedbacks and regulatory interactions forming the inflammatory-MPN model. Using mathematical modeling, we add further proof to the concept that chronic inflammation may be both a trigger of clonal evolution and an important driving force for MPN disease progression. Our findings support intervention at the earliest stage of cancer development to target the malignant clone and dampen concomitant inflammation.


Introduction
The classic chronic Philadelphia-negative myeloproliferative neoplasms (MPNs) include essential thrombocythemia (ET), polycythemia vera (PV) and primary myelofibrosis (PMF), which are acquired stem cell neoplasms [1]. Most patients live with their MPNs for decades although with a huge morbidity burden due to a high risk of thrombosis with cardiovascular complications and a massive comorbidity burden as well due to an increased propensity to develop autoimmune and chronic inflammatory diseases [2][3][4], including a 40% increased risk of second cancers [5,6]-not only after the MPN-diagnosis but also prior to the MPN-diagnosis [7]. Several years prior to the MPN-diagnosis these patients also have an increased risk of cardiovascular, autoimmune and inflammatory diseases [8,9]. Furthermore, the MPNs have an inherent risk of transformation to acute myelogenous leukemia (AML) and myelodysplastic syndrome (MDS) [10].
Another research field-mathematical modelling of cancer development-has not been applied to a similar extent within MPNs until very recently [78,79] and not in the context of investigating the concept of MPNs as "A Human Inflammation Model for Cancer Development". Mathematical modelling of cancer development has provided new insights regarding cancer initiation and progression [80][81][82][83][84][85][86][87][88][89]. In this context, mathematical modelling has a huge potential to support or disprove understanding of research data on pathogenetic factors of significance for cancer development but also in regard to providing supportive evidence for a drug to be used in cancer therapy and accordingly a novel tool in evidencebased medicine [90][91][92]. Mathematical modelling of chronic inflammation as the trigger and driver of MPNs has never been investigated. Although the concept of MPNs as "inflammatory diseases" is being increasingly recognized, additional proof of this novel concept by mathematical modelling might be of utmost importance not only for our understanding of the pathogenesis of these neoplasms, but also in regard to diagnosis and treatment. Herein, we for the first time by mathematical modelling add further proof of the concept that MPNs may be both triggered and driven by chronic inflammation. We discuss the perspectives of our findings, which might implicate intervention at the earliest stage of cancer development (ET, PV) to target the malignant clone and dampen concomitant inflammation when the tumor burden is minimal, and accordingly, the outcome of treatment is logically most favorable.

Methods
The system describes the proliferation from stem cells to mature cells including mutations of healthy stem cells to become malignant stem cells. We include regulatory interactions (e.g. niche growth effects) and inflammation coping with cell death, inflammatory cytokines, and neutrophils. In order to design an inflammatory MPN model, we build on the coupled dynamics of inflammation and cancer progression as depicted in Fig 1.

The model
Most previous studies attempting to model the role of inflammation and immune deregulation in cancer progression consider solid tumors and couple the T-cell and natural killer (NK) cell dynamics to a logistic growth of a tumor. They mainly describe quite simplified versions of the adaptive immune response without explicitly considering the underlying cancer growth dynamics [93][94][95][96][97][98]. In contrast to all these models, our model is the first which couples the principles underlying actual cell dynamics to a basal inflammatory response. This response is seen for a normal infection, where the amount of dead cells provokes the immune response and stimulates the renewal of stem cells. Despite this complex coupling, the model is kept as simple as possible still allowing the relevant quantities to be described. Thus, the goal is to describe an important coupling between MPN development and the inflammatory response at The conceptual model. Light gray boxes (symbolized x 0 , x 1 , y 0 , and y 1 ) illustrate the compartments of the basic model, and the black arrows the rates of the flows between these compartments. Here x 0 denotes the number of HSC, x 1 that of HMS, y 0 that of MPN SC, and y 1 the number of MPN MC. The light blue compartment (symbolized a) contains all dead cells and the light orange compartment (symbolized s) the inflammatory level, i.e. the immune response. Blue arrows from these represent related rates of flows. Red stipulated arrows going from the inflammatory compartment represent effects of the cytokines (or neutrophils when eliminating dead cells) modulating rates of the basic model. Two additional rates (depending on x 0 and y 0 ) appearing as red stipulated arrows represent the bone marrow niches symbiosis with the stem cells modulating the self-renewal rates. Note, stem cells leaving their respective compartments enter the corresponding mature cell-pools as multiplied by the progenitor amplification factor (A). a quantitatively conceptual level. Hence, the complicated mathematical question of model identifiability and accurate parameter estimation will be addressed elsewhere. Nevertheless, we include some model calibration and validation after presenting the model to justify and demonstrate the strength of the model. Basically, our model consists of four pools of cells; the hematopoietic stem cells (HSC), the hematopoietic mature cells (HMC), the MPN-mutated stem cells (MPN SC) and the MPN mature cells (MPN MC). The number of these cells are denoted x 0 , x 1 , y 0 , and y 1 respectively, where x refers to normal hematopoietic cells and y to MPN hematopoietic cells, while index 0 refers to stem cells and index 1 to mature cells. A single stem cell (SC) may proliferate in three ways; symmetric self-renewal (having two stem cells as offspring), asymmetric self-renewal (turning into one stem cell and one progenitor cell), and symmetric differentiation (giving rise to two progenitor cells). The progenitor cells cannot be ignored, however, we consider the progenitor cells simply as intermedia multiplication steps describing the way stem cells generate mature cells. In the model, the generations or continuum of progenitor cells will be implicitly accounted for as each stem cell will generate a number of mature cells by an amplification factor, A (= 2 k if there are k generations of progenitor cells). Feedbacks from or to progenitor cells are ignored or integrated into the other included feedbacks.
The present focus is on the ensemble of each cell type and not the individual cells; thus the governing laws will be for the pools of cells, in science denoted compartments [99]. Mathematically, the dynamics will be described by non-linear ordinary differential equations respecting conservation laws. The HSC self-renews with rate r x and the malignant MPN SC self-renews with rate r y . Furthermore, HSC may be transformed by cell division by a rate a x whereas the MPN SC does so with a rate a y . The mature cells are multiply generated, i.e. the HMC are generated with a rate a x ÁA x and the MPN MC with a rate a y ÁA y . Finally, all cell types may die; stem cells with a lower rate and mature cells with a higher rate. The turnover (or mortality) rates are d x0 , d x1 , d y0 , and d y1 for the HSC, HMC, MPN SC, and MPN MC, respectively. Except for the mutation part and the multiplication factor, this duplicates the structure of the model proposed by Dingli and Michor (they silently used A = 1) [92].
A small probability r m describes the mutation of HSC into MPN SC. In that case, r m is not the probability of a single mutation but possibly a serial sequence of mutations turning the HSC into a cancer cell capable of self-renewal, by definition an MPN SC, where a mutation is expected to be described by a Poisson process [100]. The probability for one mutation is about 10 −7 per year per cell [101]. However, not all mutations are malignant; only mutations which happen on particular locations (i.e. at specific nucleic acids) of the DNA cause MPN relevant mutations. Inflammation increases the risk of mutations, including smoking, exposure to ultraviolet light or certain chemicals [49,50, [101][102][103][104]. It is this small probability which violates a possible deterministic description with a simple mutation rate. Except for the mutation part, the model will be deterministic and continuous. In most of our work, we studied the development right after the first malignant mutation has occurred (denoted the first insult). In these cases, the simulations start with one malignant stem cell. Meanwhile, the number of all other cells are in a healthy steady state with the mutation rate put to zero. The approach is justified by the fact that including a non-zero mutation rate did not affect the outcome of the model. resulting in specific systems of ordinary differential equations as given in S1 Appendix.
Whenever cells die the debris have to be engulfed by phagocytic cells, e.g. neutrophils and macrophages while a hierarchic cascade of pro-and anti-inflammatory cytokines are released [96][97][98]110]. Following the parsimonious principle, we let the dead cells (a) up-regulate the amount of phagocytic cells (s) with rate constant r s per dead cell while they are eliminated with a rate e s . In addition, endotoxins, smoking and other environmental factors may add to the inflammatory response; thus we add such a term (characterized by the lightning symbol in Fig  1). Since MPNs develop on time-scale years and inflammatory immune processes are fast (on time-scale hours-days), we assume that the amount of phagocytic cells is balanced by the cytokines levels in a fixed ratio. Thus, the cytokine level is proportional to the phagocytic level why the inflammatory compartment (s) represents both (up to a possible proportionality constant which may be incorporated into the rate constants). Meanwhile, the amount of dead cells is down-regulated as a second order elimination process, -e a ÁaÁs, with rate constant e a . Dead cells are produced by d 0 Áx 0 +d y0 Áy 0 +d x1 Áx 1 +d y1 Áy 1 per time denoted the turnover, which is assessed by the plasma concentration of lactic dehydrogenase (LDH). It is well-known that the inflammatory level affects the mutation rate [104] and the self-renewal rates [105]. For simplicity, we take these to be proportional with the inflammatory level (of course saturation may occur) but since the level (a) settles at constant levels so does the inflammatory level (s), which may be thought of as the amount of inflammatory cytokines which have been shown to be increased in patients with MPNs and several in a step-wise manner from controls over the early cancer stages (ET, PV) to the advanced cancer stage-myelofibrosis (PMF) (S1 Appendix) [40-46]. Thus, it turns out that various specific cytokines (IL-1β, IL-1RA, IL-2R, IL-6, IL-8, IL-10, IL-12) and C-reactive protein (CRP)-a conventional biomarker of inflammation-are linearly correlated with the inflammatory level (s). These cytokines have been chosen for validation of our model since elevated levels of several of these cytokines have been associated with an inferior survival [44]. Likewise, elevated levels of CRP have been shown to be associated with  Table B in S1 Appendix including default parameter values.

Model calibration, validation, and results
The model is inspired by Dingli and Michor, and therefore the parameter values are based upon their values [92]. However, we have adjusted them to obtain more appropriate saturation levels in agreement with data (see Fig 2 and the reported values in Table C in S1 Appendix). First, the model is calibrated to the situation of no MPN cancer cells (y 0i = 0 and y 1i = 0). In this situation, we expect a stable steady state such that the number of HSC is approximately 10 4 and that of HMC is approximately 10 10 . These choices are compromises between reported values for the number of HSC [78, 86,88,89,92].
From the steady state condition we have the number of dead cells to be e a s % 10 3 . We further expect r x >d x0 + a x and d x0 (d x1 . When allowing for MPN development the healthy state becomes unstable when perturbed by the malignant stem cells. Thus, we expect r y >r x .
In the final stage the in silico patient will have vanishing hematopoietic cells and the MPN cells will approach a stable steady state with a higher amount of MPN cells than normal hematopoietic cells in the healthy steady state. This is accomplished by choosing all the c-values equal in order to keep the model as simple as possible and the number of parameters as few as possible. Likewise, the parsimonious principle suggests d y0 = d x0 , a y = a x and A y = A x .
All these attempts in calibrating the model may simultaneously be considered as validation since they performed successfully. However, the model may be validated further by predicting affected cytokine levels from the inflammatory level. As indicators of the inflammatory level, we refer to those cytokines, which are considered most important in the context of MPNs: IL-1β, IL-1RA, IL-2R, IL-6, IL-8, IL-10, IL-12 and the inflammation biomarker CRP which all turned out to be linearly correlated with the inflammatory level (s).
For the specific cytokines (C i ) tabulated in the S1 Appendix, we have 'Normal', 'PV', and 'PMF' median values (m ij , where index i specifies the cytokine and index j refers to 'Normal', 'PV' and 'PMF' states) for each. Then we find k i1 and k i2 such that m ij = k i1 s j +k i2 where s j is the value of s at year t j . Similarly, LDH values were demonstrated to be correlated and compared to the total rate of dying cells DI = dx 0 x 0 + dx 1 x 1 + dy 0 y 0 + dy 1 y 1 . The results are summarized   Having a continuous mutation rate, it will take 24 years for the disease to develop to an allele burden of 7% (e.g. ET) and after additional four years the allele burden reaches 33% (e.g. transformation of ET to PV) to become 67% (e.g. transformation of PV to post-PV myelofibrosis) at year 36 after the first stem cell mutation. Disallowing mutations in the model and initially including a single malignant stem cell and no malignant mature cells shifts the allele burden curve by one year to the left on the time axis.
Thus, the mutation of an HSC to MPN SC triggers the disease. Once an MPN stem cell is established the disease can progress without further mutations.
The baseline inflammatory load (stimulus) is arbitrarily set to 7 pg/ml per day during normal circumstances. It is an exogenous stimulation of the immune system, which leads to an inflammatory level of 3.61 pg/ml, increasing to 3.66 pg/ml in MPNs. This corresponds to a baseline of 700 dead cells (in the hematopoietic steady state) before MPN develops remarkably. A doubling of the baseline inflammatory level is directly affecting the inflammation load (cytokine level) and thereby affecting the rest of the system as dictated by the model equations. In Fig 5 is depicted that shortening the exposure time of inflammation load is associated with deceleration of disease progression.

Discussion
Chronic inflammation is characterized by persistently activated immune cells, DNA damage, tissue destruction, remodeling and fibrosis [106]. In patients with MPNs, these processes are exemplified by the advanced myelofibrosis stage [4,29], which accordingly might be considered to develop as the consequence of chronic inflammation in the bone marrow-"the inflamed bone marrow" and "the wound that won't heal" [4,29,107,108]. Herein, we for the first time use mathematical modelling to substantiate the concept that MPN progression is facilitated by chronic inflammation and that ET and PV are linked through increasing JAK2V617F allele burden [19] which is destined to happen as time increases without interference. Importantly, we were able to create the inflammation-MPN model based upon current knowledge on the interactions between inflammatory cytokines, hematopoietic stem cells and progenitors, and the bone marrow microenvironment [31-33,35-37,105]. By mathematical modelling of all these interactions, our integrated inflammation-MPN model was created. The model was validated from current data on circulating inflammatory cytokines in MPNs [40,44-46], thereby substantiating inflammation to be a highly potent stimulus for clonal evolution and cancer progression in MPNs. In the context that elevated CRP levels have been shown to be associated with shortened leukemia-free survival in myelofibrosis [42], it is of interest that our model was excellently validated by data on CRP levels in the different MPN disease stages as well.
Mathematical modelling has been used to describe the impact of chronic inflammation and immune deregulation in aging [109] and several diseases, including type 1 diabetes mellitus [110], rheumatoid arthritis [96] and colitis-associated colon cancer [111]. Based upon the known association between respiratory infections and chronic inflammation, Herald described a general model of inflammation [97]. In this model, a system of nonlinear ordinary differential equations was used to describe interactions between macrophages, inflammatory and antiinflammatory cytokines and bacteria. Though initiated by bacteria as the stimulus to trigger chronic inflammation, their study focused on chronic inflammation in the absence of pathogens as well [97]. Of note, even small changes in parameters of importance for inflammatory cytokine production and macrophage sensitivity to cytokines resulted in dramatically different model behaviors [97]. According to this model chronic inflammation is not triggered when the immune system is functioning properly. However, in patients with a dysfunction of the immune system positive feedback of the inflammatory cytokine network is prone to induce chronic inflammation. Furthermore, if the macrophage population is more sensitive to inflammatory cytokines small perturbations initiated by the inflammation stimulus will also lead to chronic inflammation [97]. In this context, it is intriguing to consider if the inherited genetic predisposition to acquire the JAK2V617F-mutation due to the haplotype 46/1 [112][113][114][115][116][117], which also confers an increased risk of (other) inflammatory diseases (e.g. Crohns' disease) [118,119] and /or acquired genetic instability due to sustained chronic inflammation (chronic inflammatory diseases or toxin exposure (e.g. smoking) might further increase the risk of developing MPN-a hypothesis originally proposed by Hermouet et al [33,35]. Importantly, the hypersensitivity of clonal MPN-cells to exogenous and endogenous growth factors and inflammatory cytokines might also more easily lead to a chronic inflammatory state-similar to the increased sensitivity of the macrophage population leading to chronic inflammation in the . Lower: Displaying the results in terms of the clinically available quantity, total blood cell count, also shows a boosted effect with increasing inflammatory load (C). The allele burden of JAK2 mutated blood cells similarly shows that increased inflammation increases disease development (D). There is a clear effect of MPN promotion with increasing inflammatory load, earlier onset, and exposure time. Lowering inflammatory load makes disease progression less rapid. Maintaining a doubling (red dotted curve) shifts the allele burden curve to the left by two years. Shortening the exposure time of inflammatory load weakens the disease progression. The inflammation has a fast impact on the total number of blood cells, which typically changes by 25% within the first year after doubling or reducing the inflammatory load by 50%. https://doi.org/10.1371/journal.pone.0183620.g005 Herald model and also implemented in the Hermouet model, implying an enhanced myelomonocytic response to cytokine stimulation [33,35].
In the Herald model and the model described by Nielsen et al in regard to type 1 diabetes mellitus, the macrophages constituted an important compartment [97,110]. The monocytemacrophage cell lineage is of major importance in the context of inflammation and cancer development. In our MPN-inflammation model bone marrow macrophages are also of utmost importance-both in regard to release of inflammatory cytokines, but also in regard to the development of myelofibrosis. Thus, in MPNs the monocyte-macrophage cell-together with the megakaryocyte (MK) cell lineage-are considered to be responsible for the development of myelofibrosis by the release of a number of growth factors and inflammatory cytokines that stimulate fibroblast proliferation [36, 120,121]. The "Herald Model" is in several aspects equivalent to our model when considering substituting "bacteria"in the "Herald Model" by any noxious inflammatory stimulus. In fact, we implement yet another cell lineage-the MKs-as the source of a continuous release of products that stimulate the vicious inflammation circle, implying ultimately the development of cancer-the MPNs. As previously outlined, our mathematical modelling of the concept of chronic inflammation in MPNs is also supported by the elegant model described by Hermouet and co-workers [33,35], in which the JAK2 46/1 haplotype was proposed as a marker of inappropriate myelomonocytic response to cytokine stimulation, leading to increased risk of inflammation, myeloid neoplasms, and impaired defense against infection [33].  [122], describing how chronic inflammatory processes might impinge on hematopoiesis, potentially fostering hematopoietic stem cell diseases, including MPNs. By inducing high proliferation of most HSCs, chronic inflammation might give rise to both exhaustion of the HSC pool and an even greater risk to accumulate genetic alterations in HSCs. Furthermore, by inflammatory stimuli from the bone marrow microenvironment these genetically altered HSCs might be rescued or "cancer cell niche" for later development of a hematological cancer [122].
The perspectives of our study are several. In the context that myelomonocytic cells (granulocytes, macrophages, monocytes) and MKs are all deeply involved in cancer development and progression [123,124], chronic inflammation is associated with premature atherosclerosis (atherothrombosis) [21-23, 29,30], in which both platelets and monocytes are highly important (monocytes a link between atherosclerosis and cancer [28]) and platelets are intimately involved in the metastatic process in cancer [124]-and likely in MPNs as well [125]-the avenue is opened for studying all these aspects by using mathematical modelling of current knowledge of the impact of chronic inflammation and immune deregulation in patients with MPNs. Ultimately, mathematical modelling may also be able to substantiate which agents to be used in MPNs in order to induce "minimal residual disease" [125][126][127][128][129] and the importance of early intervention with agents that directly target both the malignant clone (interferon-alpha2) [126-129] and the inflammatory process (JAK1-2 inhibition with e.g. ruxolitinib) [130].
In conclusion, we have for the first time applied mathematical modelling as a tool to deliver the proof of concept that chronic inflammation is closely linked to the development of the MPNs-myeloproliferative cancers which today are considered to be "chronic inflammatory diseases", in which chronic inflammation may be a driving force for clonal expansion and ultimately the development of AML [4,[29][30][31][32]39]. Studies are ongoing to elucidate the above perspectives by mathematical modelling. In this regard, mathematical modelling of resolution of inflammation may be highly important [98] and useful to support the decision-making which agents to use in the future for patients with MPNs in order to induce minimal residual disease and hopefully cure.