User:Asf45/Final Term Paper

From BIOL 300 Wiki
Jump to: navigation, search

Introduction

Multiple Myeloma Bone Disease (MMBD) is a plasma cell malignancy that disrupts the natural process of bone repair.[1] In patients with MMBD, malignant plasma cells in the bone marrow can cause bone lesions, osteoporosis, and ultimately breakage.[2][3] The progression of the disease depends on autocrine and paracrine interactions between myeloma (malignant plasma) cells, osteoclast (bone-degrading) cells, and osteoblast (bone-synthesizing) cells. Paracrine interactions become particularly important to bone degradation in the later stages of the disease as myeloma cells impede osteoblast growth to the point that is severely suppressed or even completely absent[4] and promote osteoclast growth beyond normal levels.[5] Additionally, osteoclasts have been found to promote myeloma cell growth, further accelerating the progression of the disease.[5][6] More recently, a mechanism has been proposed for the fusion of osteoclast and myeloma cells within tumors to form joint cells which allow malignant cells to participate in the degradation of bone directly. [7]


There is currently no cure for MMBD, but the illness can be controlled with treatment. [8] Because the bone marrow microenvironment is known to support tumor cell growth, it is necessary to understand the dynamics of the interactions between osteoclast, osteoblast and myeloma cells in the marrow in order to develop more effective treatments and ultimately a cure for MMBD.[5] Additionally, the inclusion of recently-discovered joint cells in models of bone cell dynamics will be essential to our understanding of possible relapse mechanisms of the disease as the consequences of such osteoclast-myeloma cell fusions is currently unclear.[3]


The model presented in this paper focuses on the importance of paracrine interactions in the progression of MMBD in a representative volume of bone marrow. In the first part of this model reconstruction, we tested the hypothesis that increased paracrine interaction intensity is a major driver of unstable bone cell destruction in patients with MMBD (Multiple Myeloma Bone Disease). This model builds on work by Komarova et al. (2003) in which autocrine interactions are identified as major drivers of bone growth dynamics in the absence of MMBD.[9] In the second part of this model reconstruction, we introduced recently discovered joint cells to the model to test the hypothesis that the shielding of myeloma cells from degradation when they are part of a joint cell alters the relapse mechanism of MMBD by changing the dynamics of relative cell population sizes over time. Finally, we extended the model by conducting a sensitivity analysis to determine to which parameters the model is most sensitive.

Replication of the model was for the most part successful, but some discrepancies between our results and those in the paper led us to believe that the list of parameter values given in the paper was incorrect. Despite these discrepancies, all replicated figures from the first part of the model were qualitatively consistent with the figures in the paper, which allowed us to confirm our first hypothesis that increased paracrine interaction intensity is a major driver of unstable bone cell destruction in patients with MMBD. Extension of the model indicated that it was sensitive to certain autocrine parameters as well, which could indicate that autocrine interaction intensity also plays an important role in unstable bone cell destruction. For the second part of the model, the discrepancies between the figures that we reproduced and the figures in the paper were larger, but despite these discrepancies it was still clear that joint cells altered bone growth dynamics when they were introduced to the model. Thus, we were able to confirm that the second part of our hypothesis proposing that the shielding of myeloma cells from degradation when they are part of a joint cell alters the relapse mechanism of MMBD by changing the dynamics of relative cell population sizes over time was correct.

Model Description

This model description outlines my reconstruction of a model describing the changes in the numbers of osteoblast, osteoclast, myeloma and joint cells in a representative volume of bone tissue in patients with Multiple Myeloma Bone Disease (MMBD) presented by Koenders & Saso (2016) in A mathematical model of cell equilibrium and joint cell formation in multiple myeloma.[3] This will include a list of the state variables and parameters, a discussion of the assumptions of the model, an explanation of the equations presented in the paper, and an explanation of the steps taken in recreating the model.

State Variables and Parameters

Table 1: State variables
Variable Description
B Number of osteoblast cells in BMU
C Number of osteoclast cells in BMU
T Number of myeloma cells in BMU
J Number of joint cells in BMU
  • Note: BMU = "basic molecular unit" consisting of a representative volume of bone marrow.
Table 2: Independent parameters
Variable Description Value
gCC Autocrine parameter representing the acuity of growth in osteoclasts due to the presence of osteoclasts 0.5
gCB Paracrine parameter representing the acuity of growth in osteoclasts due to the presence of osteoblasts -0.5
gBB Autocrine parameter representing the acuity of growth in osteoblasts due to the presence of osteoblasts 0
gBC Paracrine parameter representing the acuity of growth in osteoblasts due to the presence of osteoclasts 1
gTT Autocrine parameter representing the acuity of growth in myeloma cells due to the presence of myeloma cells 0.5
gTB Paracrine parameter representing the acuity of growth in myeloma cells due to the presence of osteoblasts 0
αB Quantitative growth rate of osteoblasts 4
αC Quantitative growth rate of osteoclasts 3
αT Quantitative growth rate of myeloma cells 0.3
αJ Quantitative growth rate of joint cells 1 • 10-3
βB Quantitative decay rate of osteoblasts 0.02
βC Quantitative decay rate of osteoclasts 0.2
βT Quantitative decay rate of myeloma cells 0.1
βJ Quantitative decay rate of joint cells 0.3
κ Dissociation constant of joint cells 0.2
SL Number of myeloma cells present in smouldering myeloma 10
Table 3: Calculated parameters
Variable Description
hCT Paracrine parameter representing the increased sensitivity of osteoclasts to proliferation in the presence of myeloma cells
hBT Paracrine parameter representing the increased sensitivity of osteoblasts to decay in the presence of myeloma cells
gTC Paracrine parameter representing the acuity of growth in myeloma cells due to the presence of osteoclasts
  • Note: parameters gXY, hXY, αX and βX are all in units of cells/day.

Assumptions

The first part of the model excluding the joint cells has only one assumption: the presence of myeloma does not influence the interaction between osteoblast and osteoclast cells. The second part of the model including the joint cells has three assumptions: 1) the increase in the number of joint cells is proportional to the amount of myeloma and osteoclast cells present; 2) when joint cells dissociate, the myeloma cells survive but the osteoclast cells do not, and 3) there is no paracrine interaction between joint cells and osteoblast or remaining osteoclast populations. Since joint cells are a relatively new discovery, their interactions with other cells are not yet fully understood, and thus this model aims to address only the effects of the formation/dissociation of joint cells on the progression of the disease.

Equations

Rate Equations

The following three equations describe the growth of osteoclast (C), osteoblast (B), and myeloma (T) cell populations in a representative volume of bone marrow referred to as a basic molecular unit (BMU).

Eq. 1:


  \frac{dC}{dt} = \alpha_{C}C^{g_{CC}}B^{g_{CB}} (1+h_{CT}T)-\beta_{C}C

This equation describes the growth of the osteoclast cell population in a BMU. The first term, αCCgCC, represents the growth of osteoclast cells as characterized by their growth rate αC and the autocrine parameter gCC. The second term BgCB accounts for the effects of paracrine interactions between osteoclast and osteoblast cells. The third term (1+hCTT) accounts for the promotion of osteoclast growth via paracrine interactions with myeloma cells. Finally, the last term -βCC accounts for osteoclast apoptosis as characterized by the natural osteoclast decay rate βC and the population size C.


Eq. 2:


  \frac{dB}{dt} = \alpha_{B}B^{g_{BB}}C^{g_{BC}} (1-h_{BT}T)-\beta_{B}B

This equation describes the growth of the osteoblast cell population in a BMU. The first term, αBBgBB, represents the growth of osteoclast cells as characterized by their growth rate αB and the autocrine parameter gBB. The second term CgBC accounts for the effects of paracrine interactions between osteoclast and osteoblast cells. The third term (1-hBTT) accounts for the suppression of osteoblast growth via paracrine interactions with myeloma cells. Finally, the last term -βBB accounts for osteoblast apoptosis as characterized by the natural osteoblast decay rate βB and the population size B.


Eq. 3:


  \frac{dT}{dt} = \alpha_{T}T^{g_{TT}}C^{g_{TC}} -\beta_{T}T

This equation describes the growth of the myeloma cell population in a BMU. The first term, αTTgTT, represents the growth of myeloma cells as characterized by their growth rate αT and the autocrine parameter gTT. The second term CgTC accounts for the effects of paracrine interactions between osteoclast and myeloma cells. Finally, the last term -βTT accounts for myeloma cell apoptosis as characterized by the natural myeloma cell decay rate βT and the population size T.


Equilibrium equations

The following three equations describe the equilibrium values for osteoclast, osteoblast, and myeloma cell populations in a BMU. These equations were obtained by setting the left-hand sides of the three rate equations above equal to zero and solving for the population numbers of the participating cells, then approximating the system with gTC=0.

Eq. 4:


  \bar{T}^{(1)} = \left ( \frac{\alpha_{T}}{\beta_{T}} \right ) ^\frac{1}{1-g_{TT}}


Eq. 5:


  \bar{C}^{(1)} = \left [ \frac{\alpha_{C} \left ( 1+h_{CT}\bar{T}^{(1)} \right )}{\beta_{C}} \right ] ^\frac{(1-g_{BB})}{ \left (g_{CC}(g_{BB}-1)-g_{CB}g_{BC}-g_{BB}+1 \right )} * 
\left [ \frac{\alpha_{B} \left ( 1-h_{BT}\bar{T}^{(1)} \right )}{\beta_{B}} \right ] ^\frac{(g_{CB})}{ \left (g_{CC}(g_{BB}-1)-g_{CB}g_{BC}-g_{BB}+1 \right )}


Eq. 6:


  \bar{B}^{(1)} = \left [ \frac{\alpha_{C} \left ( 1+h_{CT}\bar{T}^{(1)} \right )}{\beta_{C}} \right ] ^\frac{(g_{BC})}{ \left (g_{CC}(g_{BB}-1)-g_{CB}g_{BC}-g_{BB}+1 \right )} * 
\left [ \frac{\alpha_{B} \left ( 1-h_{BT}\bar{T}^{(1)} \right )}{\beta_{B}} \right ] ^\frac{(1-g_{CC})}{ \left (g_{CC}(g_{BB}-1)-g_{CB}g_{BC}-g_{BB}+1 \right )}


Rate equations including joint cells

The following three equations describe the growth of osteoclast, osteoblast, myeloma, and joint (J) cell populations in a BMU.

Eq. 7:


  \frac{dC}{dt} = \alpha_{C}(1+h_{CT}T)C^{g_{CC}}B^{g_{CB}}-\beta_{C}C-\alpha_{J}CT

This equation is the same as the rate equation for osteoclasts in the first part of the model excluding joint cells (Eq. 1) but subtracts the term αJCT. The subtraction of this term represents the subtraction of one osteoclast cell from its cohort when a joint cell is formed. There is no term added to account for the addition of a single myeloma cell to its cohort when a joint cell dissociates because it is assumed that when joint cells dissociate, only the myeloma cell survives.


Eq. 8:


  \frac{dB}{dt} = \alpha_{B}(1-h_{BT}T)C^{g_{BC}}B^{g_{BB}} -\beta_{B}B

This equation is the same as the rate equation for osteoblasts in the first part of the model excluding joint cells (Eq. 2) because it is assumed that osteoblast population growth is not affected by the addition of joint cells to the model.


Eq. 9:


  \frac{dT}{dt} = \alpha_{T}C^{g_{TC}}T^{g_{TT}}-\beta_{T}T-\alpha_{J}CT+\kappa J

This equation is the same as the rate equation for myeloma cells in the first part of the model excluding joint cells (Eq. 3) but subtracts the term αJCT and adds the term ΚJ. The subtraction of the term αJCT represents the subtraction of one myeloma cell from its cohort when a joint cell is formed, and the addition of the term ΚJ represents the addition of one myeloma cell to its cohort when a joint cell dissociates.


Eq. 10:


  \frac{dJ}{dt} = \alpha_{J}CT-\beta_{J}J

This equation is added to the model to describe the growth of the joint cell population in a BMU. The first term, αJCT, represents the formation of the joint cells as characterized by their growth rate αJ and the amount of osteoblasts C and osteoclasts B in the BMU. The second term accounts for the apoptosis of joint cells as characterized by their decay rate βJ and their population size J.

Simulation

The equations listed above were simulated using Mathematica. First, I (1) attempted to find the equilibrium solutions to the rate equations (Eqs. 1-3) using Solve, but a solution was not possible. Instead, I demonstrated that the equilibrium equations (Eqs. 4-6) were solutions to the rate equations (Eqs. 1-3) by plugging them into the rate equations along with the appropriate parameter values and showing that they reduce to zero. Next, I (2) used Plot and NDSolve with the rate equations for the three cell populations to visualize the equilibrium values for the cell populations relative to their healthy values (C0, B0, T0, and J0) under high myelomic paracrine activity (hCT=hBT=0.03) and low myelomic paracrine activity (hCT=hBT=0.01) for values of the paracrine parameter gTC ranging from 0 to 1.2. Next, I (3) used Manipulate, Plot, and NDSolve to show how varying the magnitude of an initial perturbation from equilibrium would alter cell population responses in the absence and presence of myeloma cells. I then (4) simulated cell population growth for low and high myelomic paracrine activity when joint cells were added to the model using the same procedure as in (2). Next, I (5) simulated the relapse of MMBD from a reservoir of joint cells by simulating the development of the cell populations after all myeloma cells had been eradicated using Plot and NDSolve for the initial conditions C/C0[0]=1.0, B/B0[0]=1.0, T/T0[0]=0, and J/J0[0]=1.4. Finally, I extended the model by conducting a sensitivity analysis to determine the parameters most important to the progression of MMBD.

Results

Most figures from the paper were able to be successfully replicated within a small margin of error with the exception of Fig. 7. As will be discussed below, several small discrepancies indicate that at least one parameter was reported incorrectly by the authors, leading to the differences between my replicated figures and the figures from the paper. Generally, however, the biological implications of each of my replicated figures with the exception of Fig. 7 are no different than those of the figures in the paper.

Figure 1

Figure 1 from Koenders & Saso (2016).

Replication of Figure 1 was unsuccessful because the authors did not provide enough information to replicate it. The authors describe this figure as "denoting limiting stability for various choices of the paracrine osteoclast growth stimulation by tumor cells hCT [where] the model is stable for values of the exponent gTC that are smaller than the line values." This figure most likely represents nullclines, but without knowing how the authors combined the three rate equations (Eqs. 1-3), it is impossible to construct them. Another important note here is that the authors do not indicate why they multiply the parameter hCT by ST, and without knowing this information, we cannot know how it might affect the model and subsequently the nullclines. The authors indicate earlier in their paper that the parameter SL represents the number of tumor cells in smouldering myeloma and is "of the order of 10", but do not make any reference to it at any other point in the paper. Here, we can only assume that a mistake occurred in the assembly of this figure where ST is supposed to be SL=10, because the parameter ST appears only in Figure 1.

Figure 2

Figure 2 describes the cell population sizes of osteoclast, osteoblast and tumor cells under high myelomic paracrine activity (hCTSL=hBTSL=0.3) and low myelomic paracrine activity ((hCTSL=hBTSL=0.1) relative to their healthy values (hCTSL=hBTSL=0) as a function of the paracrine parameter gTC. This parameter represents the acuity of growth in myeloma cells due to the presence of osteoclast cells. In my recreation of Figure 2, blue=osteoblast, green=osteoclast, and red=myeloma. Model parameters are αC=3, αB=4, αC=3, βC=0.2, βB=0.02, αT=0.3, βT=0.1, gCC=0.5, gBB=0, gCB=-0.5, gBC=1.0, gTT=0.5, SL=10.

Figure 2 in Koenders & Saso (2016)
Reproduction of Fig. 2 in Koenders & Saso (2016)

My figure implies the same overall result as the figure in the paper, namely that a combined increase in paracrine parameters leads to a worsening stability situation. Additionally, the sizes of the three cell populations relative to one another are the same as those in the original figure. However, my figure is slightly off from the one in the paper, because for both high and low myelomic paracrine activity, the lines extend further than they do in the actual figure from Koenders & Saso (2016). Since all equations were checked several times for errors, all parameters were exactly matched to those indicated in the paper, and altering the initial conditions did not change the figure very drastically, when I was recreating this figure I began to suspect that some parameter might be reported incorrectly in the paper. This suspicion was later confirmed when I attempted to replicate Figure 5 (see below). However, despite the fact that it is slightly off from the original figure, it does not affect the conclusions that can be made from it other than the values of gTC for which stability no longer exists.

Figure 3

Figure 3 is a time plot describing the restoration to equilibrium values of osteoclast, osteoblast and myeloma cells after a 10% perturbation from equilibrium at time t=0 when no myeloma is present (hCTSL=hBTSL=gTC=0). In my recreation of Fig. 3, blue=osteoblast, green=osteoclast, and red=myeloma. All other parameters as indicated in Fig. 2.

Fig. 3 in Koenders & Saso (2016)
Reproduction of Fig. 3 in Koenders & Saso (2016)

Successful replication of this figure indicates that, after 40-60 days, the natural bone repair process will return osteoclast and osteoblast cell population sizes to their equilibrium values after a small injury. Since myeloma is not present, the size of its population does not change, as indicated by the horizontal red line at y=1.00.

Figure 4

Figure 4 is a time plot describing the restoration to equilibrium values of osteoclast, osteoblast and myeloma cells after a 10% perturbation from equilibrium at time t=0 when myeloma is present. In my recreation of Figure 4, blue=osteoblast, green=osteoclast, and red=myeloma. All other parameters not listed here are as indicated in Fig. 2.

Figure 4 from Koenders & Saso (2015).
First attempt at recreating Figure 4 from Koenders & Saso (2016). All parameters are set for this figure as indicated by the authors in the paper; hBT=hCT=0.12, gTC=0.9.
Second attempt at recreating Figure 4 from Koenders & Saso (2016). All parameters are set for this figure as indicated by the authors in the paper except T[0] which is changed from 0.1 to 5.0.
Third attempt at recreating Figure 4 from Koenders & Saso (2016). All parameters are set for this figure as indicated by the authors in the paper except hBT and hCT which are changed from 0.012 to 0.0012 and gTC which is changed from 1.0 to 0.0.

My first attempt at replicating Fig. 4 using the parameters given in the paper (second from left figure) resulted in a figure that was clearly different from Fig. 4 in the paper. All three cell populations begin at the correct initial values, and the general behavior of the osteoblast and osteoclast populations are very similar to the figure from the paper. On my second attempt (third from left figure), I altered the initial amount of myeloma cells in from 0.1 cells to 5.0 cells. This corrects the shape of the curve for the myeloma cell population, but the final value of the myeloma cell population at the end of the time plot is still too high. Finally, on a third attempt (far right figure), I was able to achieve the correct shapes of the curves for the myeloma and osteoclast populations by changing parameter values from hBT=hCT=0.012 to hBT=hCT=0.0012 and gTC=1.0 to gTC=0.0., but the osteoblast population remains relatively unchanged over the course of the time plot instead of initially increasing then leveling back out to equilibrium.

Figure 5

Figure 5 describes the cell population sizes of osteoclast, osteoblast, tumor and joint cells under high myelomic paracrine activity (hCTSL=hBTSL=0.3) and low myelomic paracrine activity (hCTSL=hBTSL=0.1) relative to their healthy values (hCTSL=hBTSL=0) as a function of the paracrine parameter gTC. This parameter represents the acuity of growth in myeloma cells due to the presence of osteoclast cells. In my recreation of Fig. 5, blue=osteoblast, green=osteoclast, red=myeloma, and purple=joint. All other parameters as indicated in Figure 2, and new model parameters for joint cells include α<aub>J</sub>= 1 • 10-3, βJ=0.2, and κ=0.2.

Fig. 5 in Koenders & Saso (2016).
Replication of Figure 5 in Koenders & Saso (2016).

This figure demonstrates that for the joint cell model, the combined increase of paracrine parameters still leads to a worsening stability situation. Like Fig. 2, my replication of Fig. 5 is qualitatively the same as Fig. 5 in the paper, but the values of gTC for which stability no longer exists are larger.

Figure 6

Figure 6 is a time plot describing the relapse of MMBD from a reservoir of joint cells after all myeloma cells have been degraded. In my recreation of Figure 6, blue=osteoblast, green=osteoclast, red=myeloma, and purple=joint. Here hCTSL=hBTSL=0.35, gTC=0.3, and all other parameters are as indicated in Fig. 2.

Fig. 6 in Koenders & Saso (2016)
Replication of Fig. 6 in Koenders & Saso (2016)

In my replication of Fig. 6, the myeloma and osteoblast populations demonstrate very similar behavior to those in Fig. 6 from the paper, but the joint and osteoclast populations do not. The most obvious discrepancy here is that in the paper, the osteoclast and joint cell populations reach the same value at the end of the time plot whereas in my recreation the osteoclast population ends up much lower.

Figure 7

Figure 7 is a time plot describing the relapse of MMBD from a reservoir of joint cells after all myeloma cells have been degraded. In my recreation of Fig. 7, blue=osteoblast, green=osteoclast, red=myeloma, and purple=joint. Here hCTSL=hBTSL=0.35, gTC=0.34, and all other parameters are as indicated in Fig. 2.

Fig. 7 in Koenders & Saso (2016).
Replication of Fig. 7 in Koenders & Saso (2016).

My replication of Fig. 7 was the most divergent from the figure in the paper of all replicated figures. The model behavior in the replicated figure begins similarly to that of the figure from the paper, as joint cell and osteoclast populations initially drop while the osteoblast population remains relatively unchanged and the myeloma population increases. Past t=100, however, the figures become very different in that the cell populations in my figure reach and remain at equilibrium whereas the populations in the figure from the paper remain relatively constant for several hundred days then eventually become unstable around t=500.

Extension

For my model extension, I conducted a sensitivity analysis on the first part of the model to further explore the author's assertion that paracrine interaction intensity (and therefore the value of the paracrine parameters) is a major driver of bone cell destruction in MMBD. This was particularly helpful for identifying which parameters may be causing my figures to differ slightly from those in the paper. To conduct this sensitivity analysis, I calculated how much a 10% change from the parameter values given by the authors would change the final equilibrium values of each of the cell populations. The result is shown in Table 4 below.

Table 4: Sensitivity analysis on part one of the model where parameter values given by the authors were increased by 10% and the resultant change in the final equilibrium values of the osteoclast, osteoblast and myeloma cells was calculated. Parameter gBB is excluded because the authors give it the value 0 to which a 10% change cannot be applied. Largest change in final equilibrium value for each cell population is shown in red.
Parameter Percent change in C Percent change in B Percent change in T
gCC 4.08 1.82 6.50
gCB 33.8 24.3 52.7
gBC 2.30 1.54 3.22
αC 24.1 10.8 43.7
αB 7.80 4.65 12.07
αT 12.2 0.56 42.6
βC 15.7 9.83 26.0
βB 8.7 4.32 13.8
βT 6.64 0.60 24.5
gTT 20.4 0.04 72.4
hCT 3.45 1.50 5.40
hBT 1.51 1.08 1.98
gTC 2.90 0.12 10.5

It is clear from this sensitivity analysis that the paracrine interactions between osteoclasts and osteoblasts play an important role in the behavior of the cell populations in this model, as a 10% change in the parameter gCB induced the largest changes in the final equilibrium values of osteoclast and osteoblast populations (33.8% and 24.3%, respectively). It also induced the second highest change in the final equilibrium value of the myeloma cell population. This result supports the first hypothesis proposed in this study, as increased paracrine interaction intensity provided by an increase in the value of a paracrine parameter induced the largest changes in two of the three cell populations. The most interesting information to gain from this sensitivity analysis, however, is that the autocrine parameter gTT induced the largest change in the final equilibrium value of the myeloma cell population (72.4%). While our hypothesis does not state that autocrine parameters are unimportant in bone cell destruction in MMBD, the authors of the paper chose to focus on paracrine parameters alone and did not consider how varying the value of any of the autocrine parameters would affect the model. Here we can see that indeed at least one autocrine parameter, gTT, is very important to bone growth dynamics in MMBD as a small change in its value has been shown to induce a large change in the amount of cancerous cells in the BMU.

Discussion

Hypothesis

Koenders & Saso (2015) hypothesized that increased paracrine interaction intensity is a major driver of unstable bone cell destruction in patients with Multiple Myeloma Bone Disease (MMBD). The results presented here indicate that this is true, but that the model seems to be very sensitive to certain autocrine interactions as well. By replicating the equations for the four different cell populations and using them to simulate restoration of the cell populations to their equilibrium values after a perturbation as well as relapse mechanisms of the disease, it is clear that the paracrine parameters gTC, hCT and hBT are shown to drive the stability of the model as well as the duration of bone repair processes and relapse time. Through the model extension presented here, however, it has been shown that one autocrine parameter corresponding to myeloma cells (gTT) is a major driver of the model's stability as well.

Discrepancies

While Fig. 3 and (for the most part) Fig. 2 were successfully replicated using the equations and parameter values presented by Koenders & Saso (2015), there were discrepancies ranging from minor to substantial all other figures I attempted to replicate. The first discrepancy between the model presented in the paper and my recreation appeared in Fig. 2. The ratios of the three cell populations to one another under both high and low myelomic paracrine activity in my recreation is the same as the figure in the paper, but the values of gTC at which the system loses stability are slightly larger; in the paper, gTC≈0.35 for high and gTC≈1.10 for low myelomic paracrine activity while gTC≈0.70 for high and gTC≈1.40 for low myelomic paracrine activity in my recreation. While these values are substantially different from one another, the ratios of the cell populations to one another before they lose stability is the more important piece of information to obtain from this figure for building a qualitative picture of what is actually going on in the system. Thus I agree with the authors on their qualitative interpretation of the figure which argues that the combined increase of paracrine activity associated with all three parameters leads to a worsening stability situation, and that MMBD can manifest itself in multiple ways depending on the paracrine mechanism that dominates.

The second discrepancy I found between the model presented in the paper and my recreation appeared in Fig. 4. While all three cell populations began at the correct initial values, my first recreation was quite clearly different from Fig. 4 in the paper in that the myeloma cell population continued to increase beyond T/T0=1.10 instead of initially increasing to a maximum around T/T0=1.05 and then decreasing to a final equilibrium value around T/T0=1.02 as in the paper. Additionally, the values of the osteoclast and osteoblast populations begin to diverge rather than converge toward the end of the time plot. In an attempt to discover the source of this discrepancy, I created a Manipulate function in Mathematica to vary the parameter values and initial conditions for the time plot. By first varying the initial conditions for the differential equations of the cell populations, I discovered that altering the initial value of the myeloma cells from 0.1 to 5.0 resulted in a figure that appeared more similar to the one in the paper in that the myeloma cell population settled to a final equilibrium value rather than increasing outside of the margins of the figure (see figure second from the right in Results, Figure 4). While this alteration did fix the problem of the values of the osteoblast and osteoclast populations diverging rather than converging toward the end of the time plot that appeared in my first attempt at recreating the figure, this figure still differs from Fig. 4 in the paper in that the myeloma cell population initially increases to a maximum value around T/T0=1.06 rather than 1.05 and settles to a final equilibrium value around T/T0=1.04 rather than 1.02. By varying the values of all parameters in the model using my Manipulate function, I then discovered that altering the values of the paracrine parameters hCT and hBT from 0.012 to 0.0012 and gBC from 1.0 to 0.0 resulted in a figure in which the behavior of the myeloma and osteoclast cell populations were very nearly identical to that of the same cell populations in the figure, but in which the osteoblast cell population stayed constant throughout the time plot rather than initially decreasing to around B/B0=1.01 then decreasing back to a final equilibrium value around B/B0=1.00. Since altering the paracrine parameter gBC had such a dramatic effect on the time plot, I decided to investigate the parameter values given in Komarova et al. (2003), the paper which the authors cite as a basis for their own paper, and found that while all other parameter values between the two papers were consistent, Komarova et al. (2003) gave the values gBC=-0.5 and gCB=1.0, while the authors gave the values gBC=1.0 and gCB=-0.5. Switching these two parameter values in the model given here produces a non-real result, so the discrepancies in Fig. 4 or any other figure cannot be attributed to this one issue, but it could be that instead there is some combination of parameters that, like gBC, are set to some value other than those reported by the authors to produce the figures seen in the paper.

In replicating Figures 5, 6 and 7, similar discrepancies were found where the minimum or maximum values for various cell populations were slightly off, or the shape of a single curve did not match the paper exactly. The analysis I conducted in my recreation of Fig. 4 leads me to believe that there are at least two parameter values that are reported incorrectly, and it would be incredibly difficult for me to know exactly which parameters these might be without contacting the authors directly. Because the general form of the curves in the time plots and the relative sizes of the four cell populations in the figures I recreated are consistent with those in the paper, however, it is reasonable to conclude that the authors' overall interpretations of each figure are justified.

Part of what made it difficult to trace the sources of these discrepancies was that the response variable for every figure in the paper, with the exception of Fig. 1, was plotted as a ratio of the sizes of the cell populations in a patient with MMBD to their "healthy" values for which gTC, hCT and hBT=0. Because the authors do not explicitly report initial conditions for any simulation presented in their paper, it was difficult to determine what initial conditions should be passed to the function NDSolve when replicating the figures. The authors do provide that the ratio of osteoblasts to osteoclasts at any one time is approximately 100 to 1; however, like the ratios given as the response variables in each figure, this information does not help to determine the initial value of the osteoblast, osteoclast, myeloma and joint cell populations at time=0 for simulation purposes.

Limitations

This model is limited in that it most likely oversimplifies the complex interactions occurring between osteoclast, osteoblast, myeloma and joint cells in the bone marrow in a patient with MMBD. For example, the first part of the model excluding the joint cells assumes that the presence of myeloma does not influence the interaction between osteoblast and osteoclast cells. While this assusmption may be useful for modeling purposes, it may not actually be true. The assumptions we make for the second part of the model make it even more limiting than the first part; for example, we assume that the increase in the number of joint cells is proportional to the amount of myeloma and osteoclast cells present, because joint cells are composed of myeloma and osteoclast cells. This may not necessarily be true, as the number of joint cells could be dependent only on myeloma or osteoclast cell population size, or could be proportional to the myeloma and osteoclast populations in a ratio other than 1:1. We also make the two assumptions here that when joint cells dissociate, the myeloma cells survive but the osteoclast cells do not, and that there are no paracrine interactions between joint cells and osteoblast or remaining osteoclast populations. These three assumptions make the second part of the model very limited as to its predictive power in situations other than estimating relapse mechanisms. For example, since we do not have information on the nature of the autocrine or paracrine interactions affecting joint cells, we cannot see how the model changes with varying growth parameters relating to joint cells as we can with growth parameters relating to osteoclast, osteoblast or myeloma cells in the first part of the model. Thus, since joint cells are a relatively new discovery[3], their interactions with other cells are not yet fully understood, and thus this model is only able to address the effects of the presence of joint cells on the progression of the disease and its relapse.

References

  1. Silbermann, R., & Roodman, G. D. (2010). Clinical presentation of myeloma bone disease. In Myeloma bone disease (pp. 1-13). Humana Press.
  2. Rajkumar, S. V., & Kyle, R. A. (2005, October). Multiple myeloma: diagnosis and treatment. In Mayo Clinic Proceedings (Vol. 80, No. 10, pp. 1371-1382). Elsevier.
  3. 3.0 3.1 3.2 3.3 Koenders, M. A., and R. Saso. "A mathematical model of cell equilibrium and joint cell formation in multiple myeloma." Journal of theoretical biology 390 (2016): 73-79.
  4. Esteve, Flavia R., and G. David Roodman. "Pathophysiology of myeloma bone disease." Best Practice & Research Clinical Haematology 20, no. 4 (2007): 613-624.
  5. 5.0 5.1 5.2 D'Souza, S., Kurihara, N., Shiozawa, Y., Joseph, J., Taichman, R., Galson, D. L., & Roodman, G. D. (2012). Annexin II interactions with the annexin II receptor enhance multiple myeloma cell adhesion and growth in the bone marrow microenvironment. Blood, 119(8), 1888-1896.
  6. Cui, Q., Shibata, H., Oda, A., Amou, H., Nakano, A., Yata, K., ... & Harada, T. (2011). Targeting myeloma–osteoclast interaction with Vγ9Vδ2 T cells. International journal of hematology, 94(1), 63-70.
  7. Andersen, T. L., Boissy, P., Sondergaard, T. E., Kupisiewicz, K., Plesner, T., Rasmussen, T., ... & Delaissé, J. M. (2007). Osteoclast nuclei of myeloma patients show chromosome translocations specific for the myeloma cell clone: a new type of cancer–host partnership?. The Journal of pathology, 211(1), 10-17.
  8. Ocio, Enrique M., Paul G. Richardson, S. Vincent Rajkumar, Antonio Palumbo, Maria Victoria Mateos, R. Orlowski, Shaji Kumar et al. "New drugs and novel mechanisms of action in multiple myeloma in 2013: a report from the International Myeloma Working Group (IMWG)." Leukemia 28, no. 3 (2014): 525-542.
  9. Komarova, S. V., Smith, R. J., Dixon, S. J., Sims, S. M., & Wahl, L. M. (2003). Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling. Bone, 33(2), 206-215.

Mathematica Notebook

Mathematica Notebook