A virtual approach to evaluate therapies for management of multiple myeloma induced bone disease.

Multiple myeloma bone disease is devastating for patients and a major cause of morbidity. The disease leads to bone destruction by inhibiting osteoblast activity while stimulating osteoclast activity. Recent advances in multiple myeloma research have improved our understanding of the pathogenesis of multiple myeloma-induced bone disease and suggest several potential therapeutic strategies. However, the effectiveness of some potential therapeutic strategies still requires further investigation and optimization. In this paper, a recently developed mathematical model is extended to mimic and then evaluate three therapies of the disease, namely: bisphosphonates, bortezomib and TGF-β inhibition. The model suggests that bisphosphonates and bortezomib treatments not only inhibit bone destruction, but also reduce the viability of myeloma cells. This contributes to the current debate as to whether bisphosphonate therapy has an anti-tumour effect. On the other hand, the analyses indicate that treatments designed to inhibit TGF-β do not reduce bone destruction, although it appears that they might reduce the viability of myeloma cells, which again contributes to the current controversy regarding the efficacy of TGF-β inhibition in multiple myeloma-induced bone disease.

Introduction

Multiple myeloma (MM), a haematological malignancy developed in the bone marrow, is the most common cancer involving bone and the second most prevalent cancer involving blood cells 1, 2. Bone disease is a major complication of MM and is a significant cause of morbidity in MM patients. Up to 60% of MM patients suffer a fracture during the course of the disease, and MM induced bone destruction rarely heals. Recent research into MM bone disease has revealed that the interaction between MM cells and the bone microenvironment plays an important role in the development of the condition, and a ‘vicious cycle’ of myeloma development and bone destruction is established 2, 3, 4.

Currently, several therapies are proposed to treat MM‐induced bone disease including bisphosphonates, bortezomib and TGF‐β inhibition 2, 5, 6. Bisphosphonate treatments target high turnover skeletal sites, binding to the mineralized bone matrix within these sites 7, 8, 9. After their internalization by bone‐resorbing osteoclasts, bisphosphonates inhibit further osteoclast activity and bone resorption by suppressing the differentiation of osteoclast precursors into mature osteoclasts, promoting osteoclast apoptosis and disrupting osteoclast function 8, 9. Although bisphosphonates are already a first‐line treatment for MM‐induced bone disease 7, 10, further investigation is required to determine whether bisphosphonates have an anti‐tumour effect. Several preclinical and clinical studies suggest that bisphosphonates may either have a direct or indirect anti‐tumour effect 11, 12, 13, 14, 15, 16, 17, 18. However, other studies provide contradictory evidence and suggest that bisphosphonates do not improve patient mortality 19, 20, 21.

Suppression of bone‐forming osteoblasts, which can occur from the blockade of the differentiation of osteoblast precursors into mature osteoblasts, promotes the growth of myeloma cells as well as bone destruction through supporting the production of anti‐apoptotic factors and growth factors for MM cells 2, 22. Thus, it is suggested that stimulation of osteoblast differentiation may reduce tumour burden and bone destruction in MM patients 2, 23. Bortezomib, a boron‐containing compound with the potential of enhancing osteoblast proliferation and bone formation in MM patients, has therefore been proposed as a potential therapeutic for MM‐induced bone disease 24, 25.

TGF‐β is reported to contribute to the progression of MM‐induced bone disease 5. It is released with bone resorption and stimulates the production of osteoblast progenitors but inhibits the differentiation of mature osteoblasts. It therefore suppresses bone formation and indirectly promotes the progression of MM cells (immature osteoblast cells facilitate the growth and survival of MM cells, while mature cells enhance apoptosis of MM cells). Thus, the suppression of TGF‐β is proposed as a new approach to treat MM‐induced bone disease 5; however, some controversies still exist, and further investigation is required to confirm its potential.

In the paper, a mathematical model we described previously 3 was extended to simulate these three different strategies and determine their efficacies in MM (Tables 1 and 2).

Table 1

Definitions of the π functions. See Table 2 for definitions of RANKL, OPG, TGF‐β, PTH, IL6, SLRPs, VLA4 and VCAM1.

πact,OB_uTGFβ=TGFβK_D1,TGFβ+TGFβ
πrep,OB_pTGFβ=11+TGFβ/K_D2,TGFβ
πact,OC_aTGFβ=TGFβK_D3,TGFβ+TGFβ
πact,OC_pRANKL=RANKLK_D,RANKL+RANKL
πact,MMIL6=IL6IL6+K_D,IL6,MM,act
πact,MMVCAM1=VCAM1VCAM1+K_D,VCAM1,MM,act
πrep,OB_pVCAM1=11+VCAM1/K_D,VCAM1,OB_p,rep
πact,OB_aVCAM1=VCAM1VCAM1+K_D,VCAM1,OB_a,act
πrep,MMSLRPs=11+SLRPs/K_D,SLRPs,MM,rep
πact,RANKLPTH=PTHK_D1,PTH+PTH
πrep,OPGPTH=11+PTH/K_D2,PTH
πact,RANKLIL6=IL6IL6+K_D,IL6,RANKL,act
πact,IL6VLA4=VLA4VLA4+K_D,VLA4,IL6,act
πact,IL6TGFβ=TGFβTGFβ+K_D,TGFβ,IL6,act

Table 2

Definitions of the concentrations of RANKL, OPG, TGF‐β, PTH, IL6, SLRPs, VLA4 and VCAM1. (RANKL, receptor activator of nuclear factor kappa‐B ligand; OPG, osteoprotegerin; TGF‐β, transforming growth factor‐beta; PTH, parathyroid hormone; IL‐6, interleukin‐6; SLRP, small leucine‐rich proteoglycan; VLA‐4, very late antigen‐4; VCAM‐1, vascular cell adhesion molecule 1).

RANKL	P_RANKL,d+β_RANKL⋅OB_p1+K_A,OPG⋅OPG+K_A,RANK⋅RANK⋅β_RANKLR^RANKL⋅πact,RANKLIL6⋅πact,RANKLPTH+D_RANKL
OPG	P_OPG,d+β_OPG⋅OB_a⋅πrep,OPGPTHβ_OPG⋅OB_a⋅πrep,OPGPTHOPG_max+D_OPG+D_OPG,MM⋅MM
TGFβ	α·K_res·OC_a+S_TGFβD˜_TGFβ
PTH	β_PTH+P_PTH,dtD˜_PTH
IL6	P_IL6,d+β_IL6⋅OB_u⋅πact,IL6TGF⋅πact,IL6VLA4β_IL6⋅OB_u⋅πact,IL6TGF⋅πact,IL6VLA4IL6_max+D_IL6
SLRPs	β_SLRPs⋅OB_a+P_SLRPs,dtβ_SLRPs⋅OB_aSLRPs_max+D˜_SLRPs
VLA4	P_VLA4,d+β_VLA4·MM1+K_A,VCAM1·VCAM1_tot⋅β_VLA4R^VLA4+D_VLA4
VCAM1	VCAM_tot1+K_A,VCAM1⋅VLA4

Model Development

The mathematical model which simulates the pathogenesis of MM‐induced bone disease consists of the following key equations 3: where OB, OB, OC, MM and BV are the populations of osteoblast precursors, active osteoblasts, active osteoclasts, active MM cells and bone volume respectively, and is the variation of OB with time, for example. Eqs. (1), (2), (3), (4), (5) describe the temporal variations in concentrations of OB, OB, OC, MM and BV respectively. ‘Hill functions’ are used to represent the cellular interaction via the single ligand to receptor binding, and are denoted by π functions. The definitions of the π functions in the model equations above, and the definitions and values of the model parameters are lengthy and described in detail in the work of 3 (Open Access), but are summarized here for convenience. Table 3 contains the definitions and values of the model parameters. Any unknown parameters (i.e. those parameters where experimental data are unavailable or those which have no direct biological meaning) are calculated via a genetic algorithm (GA) as indicated in Table 3, and described in detail in 3. Briefly, because a parameter may be directly or indirectly related with one or more of the initial values of cell concentrations (listed in Table 4), e.g. and involve experimental data of the initial concentration of OBp in Table 4, these initial values are set as targets for the parameter fitting. The calculation of the model parameters is then achieved by trying different values in a domain and then selecting those that provide the best fit with the corresponding experimental data. Based on these values, the remaining unknown model parameters are then calculated according to relevant experimental data through the genetic algorithm. Thus the GA approach effectively considers all possible combinations of the unknown parameters and predicts the optimal values, as described in 3. This takes many hours on a powerful PC, potentially considering billions of combinations in its search for the optimum set. The simulation was carried using the Matlab computational software package (v7.7.0, Mathworks, Natick, USA).

Table 3

Definitions and values of model parameters used in the model of MM‐induced bone disease. (GA = genetic algorithm).

Parameters	Description	Value
D_OB_u	Differentiation rate of osteoblast progenitors	3.24e + 2/day (estimated)
D_OB_p	Differentiation rate of osteoblast precursors	3.67e − 1/day (estimated)
A_OB_a	Rate of elimination of active osteoblasts	3.00e − 1/day 32
D_OC_p	Differentiation rate of osteoclast precursors	1.73e − 1/day (estimated)
A_OC_a	Rate of elimination of active osteoclasts	1.20/day 32
K D1,TGFβ	Activation coefficient related to growth factors binding on OBu	4.28e − 4 pM (calculation by GA)
K D2,TGFβ	Repression coefficient related to growth factors binding on OBp	2.19e − 4 pM (estimated)
K D3,TGFβ	Activation coefficient related to growth factors binding on OCa	4.28e − 4 pM 32
K D1,PTH	Activation coefficient for RANKL production related to PTH binding	2.09e + 1 pM (calculation by GA)
K D2,PTH	Repression coefficient for OPG production related to PTH binding	2.21e − 1 pM 32
K D,TGFβ,IL6,act	Half‐maximal concentration of TGF‐β on promoting the production of IL‐6	1.2e − 4 pM (calculation by GA)
K D,IL6,RANKL,act	Half‐maximal concentration of IL6 on promoting the production of RANKL	0.2 pM (calculation by GA)
K D,RANKL	Activation coefficient related to RANKL binding to RANK	4.12e + 1 pM (estimated)
α	TGF‐β content stored in bone matrix	1.00 pM/% 32
D˜_TGFβ	Rate of degradation of TGF‐β	2.00e + 2/day 33
βPTH	Rate of synthesis of systemic PTH	9.74e + 2 pM/day 34
D˜_PTH	Rate of degradation of PTH	3.84e + 2/day 34
βIL6	Rate of synthesis of IL6 per cell	1.20e + 7/day 35, 36
DIL6	The degradation rate of IL6	4.99e + 1/day 37
IL6max	The maximum concentration of IL‐6	8.04e − 1 pM 38
βOPG	Minimum rate of production of OPG per active osteoblast	5.02e + 6/day (estimated)
D˜_OPG	Rate of degradation of OPG	4.16/day 39
OPGmax	Maximum possible OPG concentration	7.98e + 2 pM 40
βRANKL	Production rate of RANKL per cell	8.25e + 5/day (estimated)
D˜_RANKL	Rate of degradation of RANKL	4.16/day 41
R RANKL	Maximum number of RANKL on the surface of each osteoblastic precursor	3.00e + 6 32
RANK	Fixed concentration of RANK	1.28e + 1 pM 32
K A,OPG	Association rate constant for RANKL binding to OPG.	5.68e − 2/pM 42
K A,RANK	Association rate constant for RANKL binding to RANK.	7.19e − 2/pM 42
K res	Relative rate of bone resorption (normalized with respect to normal bone resorption)	2.00e + 2%/(pM day) 43
K form	Relative rate of bone formation (normalized with respect to normal bone resorption)	3.32e + 1%/(pM day) (calculation by GA)
D MM	MM proliferation controlled by IL‐6 and BMSC‐MM adhesion	5.50e − 2/day (estimated)
A MM	Rate of elimination of active MM cells	2.00e − 3/day 44
MM max	Maximum possible MM cell concentration	1.98 pM 45
K D,VCAM1,MM,act	Half‐maximal concentration of VLA‐4 on promoting the MM cells production	1.5667e − 4/pM (calculation by GA)
K D,VLA4,IL6,act	Half‐maximal concentration of VLA − 4VLA − 4 on promoting the IL‐6 production	1.88e + 4/pM (calculation by GA)
K D,IL6,MM,act	Half‐maximal concentration of IL‐6 on promoting the MM cells production	1.2151e − 5 pM (calculation by GA)
K D,SLRPs,MM,rep	Half‐maximal concentration of SLRPs on promoting the MM cells production	1.306e + 9 pM (calculation by GA)
K_D,VCAM1,OB_p,rep	Half‐maximal concentration of VCAM‐1 on repressing the differentiation of OB p	1.4e − 1 pM (calculation by GA)
K_D,VCAM1,OB_a,act	Half‐maximal concentration of VCAM‐1 on promoting the apoptosis of OBa	2.2e − 1 pM (calculation by GA)
βVLA4	Rate of synthesis of VLA‐4 per cell	2.04e + 6/day (estimated)
D˜_VLA4	Rate of degradation of VLA‐4	1.5/day (estimated)
R VLA4	Maximum number of VLA‐4 expressed on the surface of MM cells	5.6e + 4 46
VCAM1tot	Total concentration of VCAM‐1	1.92 pM 46
KA,VCAM1	The association rate for VLA‐4 binding to VCAM‐1	8.3e − 2/pM 47
DOPG,MM	The degradation rate of OPG by MM cells	4.16/(pM day) (estimated)

Table 4

The initial values of cell concentrations in the model.

Variables	Description	Values
OB u	Uncommitted osteoblastic progenitors	3.27e − 6 pM
OB p	Osteoblast precursors	7.67e − 4 pM
OB a	Active osteoblasts	6.39e − 4 pM
OC p	Osteoclastic precursors	1.28e − 3 pM
OC a	Active osteoclasts	1.07e − 4 pM
MM	Active MM cells	3.26e − 1 pM

Notes: MM cell concentration is at day 51; other cell concentrations are at day 1.

Modelling bisphosphonates treatment

Bisphosphonate treatments inhibit bone resorption by suppressing the differentiation of mature osteoclasts as well as promoting the apoptosis of osteoclasts. Eq. (3) describes the variation of osteoclasts with time for patients with MM‐induced bone disease. In order to investigate the efficacy of biphosphonate treatments against the disease, a parameter F.Bi, representing the degree that the bisphosphonates inhibit bone resorption, is added in Eq. (3), and the new equation is as follows:

The value of parameter F.Bi is in the range of [0, 1] and is negatively correlated to the concentration of bisphosphonate during the treatment, thus a small value of F.Bi corresponds to a large dosage of bisphosphonate, which would produce a corresponding decrease in the differentiation rate of mature osteoclasts and thus osteoclast apoptosis is stimulated. For example, when F.Bi is set as 0.7, the differentiation rate of active osteoclasts is decreased to 70% (0.7), while the apoptosis of osteoclasts increases by 30% (1 − F.Bi).

Modelling bortezomib treatment

Bortezomib stimulates osteoblast proliferation and bone formation in MM patients, which can potentially inhibit the growth of myeloma cells as well as bone destruction. Eq. (2) represents the temporal variation of osteoblasts under the condition of MM‐induced bone disease. In order to simulate bortezomib treatment, a parameter F.Bo, which represents the degree by which osteoblast differentiation is promoted, is introduced to extend Eq. (2), and the new equation is as follows:

The value of parameter F.Bo is in the range of (1, + ∞), and is positively related to the dosage of bortezomib during the treatment. For example, when F.Bo is set to 2.0, osteoblast activity is increased two‐fold.

Modelling TGF‐β inhibition treatment

TGF‐β stimulates the production of osteoblast progenitors while inhibiting the differentiation of mature osteoblasts as shown in Eqs. (1) and (2), and thus the inhibition of TGF‐β indirectly suppresses the progression of MM cells, because immature osteoblast cells facilitate the growth and survival of MM cells, while mature cells enhance apoptosis of MM cells. In addition, TGF‐β can also promote the apoptosis of active osteoclasts as shown in Eq. (3). In the model of 3, the concentration of TGF‐β is defined as follows:

The definitions and values of parameters in Eq. (8) are included in Table 3. In order to examine the potential of TGF‐β inhibition treatment against MM‐induced bone disease, a parameter F.T, which describes the degree by which TGF‐β is suppressed, is added into Eq. (8), so that the concentration of TGF‐β is updated to: where the value of parameter F.T is in the range of [0, 1], and is negatively related to the concentration of TGF‐β during the treatment; for example, an F.T value of 0.9 represents a reduction in TGF‐β concentration to 90% of its normal value. As a result, model Eqs. (1), (2), (3), which contain TGF‐β, are all updated as well.

Simulation and Analysis

In all the following simulations, MM cell invasion is assumed to occur at day 51 with the different interventions applied at day 301, once the MM and bone cell populations (OB, OB and OC) have stabilized again to their final steady state value (to 578%, 293%, 199% and 249% of their values on day 50). Because the actual relationships between parameters (F.Bi, F.Bo and F.T) and equivalent drug dosages are currently not known, sample values of the parameters are investigated initially, together with further simulations to determine the sensitivity of the results to those values.

Simulation of bisphosphonates treatment

Figures 1 to 3 demonstrate how a bisphosphonate therapy, with F.Bi = 0.7, might influence cell concentrations and bone volume. Figure 1 shows a rapid increase in the population of MM cells after their initial appearance at day 51. Treatment at day 301 then leads to a reduction in peak concentration of 16% by day 450 and a continued decrease until a stable value is achieved at day 1743 (not shown) of 4.43 times the original value (at day 51). Bone cell concentrations similarly increase with the presence of the MM cells, but quickly return to a new stable state of typically 110% of their normal values (i.e. before the invasion of the myeloma cells).

Figure 1

The variation of normalized cell concentrations with respect to their initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with bisphosphonate therapy.

Figure 3

The variation of normalized bone volume with respect to its initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with bisphosphonate therapy.

The ratio of bone cells, OBa:OCa, is predicted to decrease quickly to 80% of its initial value after invasion of the MM cells (Figure 2), leading to a linear decrease in bone volume until application of the bisphosphonate therapy at day 301 (Figure 3). At this point, the OBa:OCa ratio peaks but quickly returns to 98% of its initial value, resulting in a significant slowdown in bone destruction as shown in Figure 3 after day 301.

Figure 2

The variation of normalized ratio of OBa:OCa with respect to its initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with bisphosphonate therapy.

Figures 4 to 6 show how MM cell concentration, bone volume and OBa:OCa ratio vary for sample F.Bi bisphosphonate inhibition values (of 0.7, 0.5 and 0.3) for the same treatment strategy. In relation to their peak values at day 300, MM cell population decreases to 86.8%, 85.2% and 84% (Figure 4), and OBa:OCa ratio increases to 123%, 134% and 144% when F.Bi is set to 0.7, 0.5 and 0.3 respectively (Figure 5). As illustrated in Figure 6, when F.Bi is set as 0.7, bone destruction continues, although its rate is decreased markedly, because of the increased OBa:OCa ratio. However when F.Bi is set to 0.5 or 0.3, bone destruction is halted and bone volume begins to increase again. Thus, the simulation results suggest that a smaller value of F.Bi produces a more significant inhibition of MM cell viability and bone destruction.

Figure 4

The variation of normalized MM cell concentration with respect to the value at day 300 after use of the bisphosphonates therapy with different values of F.Bi.

Figure 6

The variation of normalized bone volume with respect to its initial value after use of the bisphosphonate therapy with different values of F.Bi.

Figure 5

The variation of normalized ratio of OBa:OCa with respect to the value at day 300 after use of the bisphosphonate therapy with different values of F.Bi.

Simulation of bortezomib treatment

Simulation results for the bortezomib therapy, applied at day 301 with F.Bo set to 2.2, are shown in Figures 7 to 9, and again present the variations in cell concentrations, bone volume and the ratio of OBa:OCa. The bortezomib causes a decrease in the population of MM cells (Figure 7), with concentrations of OBp, OBa and OCa also decreasing and approaching new equilibrium points by day 450. For OBa and OCa these levels are near their initial values before the invasion of MM cells, but OBp values are reduced by 51%. Figure 8 shows that further MM‐induced bone loss is prevented after a short period of fluctuation through the intervention with bortezomib, while the OBa:OCa ratio (shown in Figure 9) again undergoes a short period of fluctuation and then returns to a level similar to its original value without the myeloma cells. This explains the termination or inhibition of MM‐induced bone loss because of bortezomib.

Figure 7

The variation of normalized cell concentrations with respect to their initial values during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with bortezomib therapy.

Figure 9

The variation of normalized ratio of OBa:OCa with respect to its initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with bortezomib therapy.

Figure 8

The variation of normalized bone volume with respect to its initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with bortezomib therapy.

Sensitivity of the simulations to the value of F.Bo is explored in Figures 10 to 12, which show the variations in the output data with F.Bo values of 2.0, 2.2 and 2.4. Overall the simulations are not particularly sensitive to this level of variation. MM cell population decreases to 86.4%, 86.2% and 86.0%, bone volume decreases to 97.4%, 97.5% and 97.6% and OBa:OCa ratio increases to 122.1%, 124.5% and 126.7%, respectively. In Figure 11, when F.Bo equals 2.0, MM‐induced bone loss continues although its rate is greatly reduced, because of the increased OBa:OCa ratio (Figure 12). When F.Bo is 2.2 and greater, a near zero or positive bone balance is achieved after the bortezomib therapy. The results suggest that the degree of reduction in MM cell viability and mitigation of bone destruction are both positively related to the value of F.Bo.

Figure 10

The variation of normalized MM cell concentration with respect to the value at day 300 after use of the bortezomib therapy with different values of F.Bo.

Figure 12

The variation of normalized ratio of OBa:OCa with respect to the value at day 300 after use of the bortezomib therapy with different values of F.Bo.

Figure 11

The variation of normalized bone volume with respect to its initial value after use of the bortezomib therapy with different values of F.Bo.

Simulation of TGF‐β inhibition treatment

The variations in cell concentrations, bone volume and the ratio of OBa:OCa after the intervention of TGF‐β therapy (from day 301 with F.T set to 0.7) are presented in Figures 13 to 15. TGF‐β suppression leads to a decline in MM cell population and bone cell concentrations, with the reduction of 13.0% in MM cell numbers (at day 450) suggesting that the tumour burden can be reduced through the inhibitionof TGF‐β. However, the MM‐induced bone destruction actually increases after the TGF‐β therapy (Figure 14). This increase in bone loss can be explained by the 18.7% decrease in OBa:OCa ratio, compared to the value at day 300, caused by the TGF‐β therapy (as shown in Figure 15).

Figure 13

The variation of normalized cell concentrations with respect to their initial values during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with TGF‐β therapy.

Figure 15

The variation of normalized ratio of OBa:OCa with respect to its initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with TGF‐β therapy.

Figure 14

The variation of normalized bone volume with respect to its initial value during different periods. Days 1 to 50: normal period; days 51 to 300: invasion of MM cells; and from day 301: intervention with TGF‐β therapy.

Finally, Figures 16 and 17 show the variations in the output data with different values of F.T (of 0.7, 0.8 and 0.9). As a result, the population of MM cells decreases to 87.0%, 87.4% and 87.8%, and the OBa:OCa ratio decreases to 81.3%, 89.6% and 97.41% of the normal value, respectively. The decrease in the ratio of active osteoblasts to osteoclasts (observed with all applications of TGF‐β) leads to continued bone loss after treatment, with the lowest levels of TGF‐β leading to the greatest loss of bone volume (Figure 18).

Figure 16

The variation of normalized MM cell concentration with respect to the value at day 300 after use of the TGF‐β therapy with different values of F.T.

Figure 17

The variation of normalized ratio of OBa:OCa with respect to the value at day 300 after use of the TGF‐β therapy with different values of F.T.

Figure 18

The variation of normalized bone volume with respect to its initial value after use of the TGF‐β therapy with different values of F.T.

Discussion

In this paper, a mathematical model first described in 3 was extended to simulate and evaluate three different therapeutic approaches to manage MM‐induced bone disease. Full details of the basic model and its validation are described in length in that publication, and are not repeated here. The therapies investigated are: bisphosphonates, bortezomib and TGF‐β inhibition, and their effects on MM and bone cell populations and bone volume are considered.

Bisphosphonates are used widely in the management of MM‐induced bone disease, and are able to inhibit osteoclast activity and bone resorption. However, the degree to which it affects MM cell viability and has an anti‐tumour effect is not clear. The model simulation suggests that bisphosphonate therapy can not only suppress bone loss, but also reduce MM cell population. This is confirmed by published data that reporting that bisphosphonates suppress MM‐induced bone destruction 8, 9. It is should be noted that direct anti‐tumour effects from the bisphosphonate are not considered in the model; thus, the decreased tumour burden is due solely to the inhibited osteoclast activity, indicating that bisphosphonate therapy has an indirect anti‐tumour effect. This finding agrees with experimental observations that a decrease in osteoclast activity inhibits proliferation of MM cells 6, 26. The underlying mechanism for the indirect effect of bisphosphonate treatment lies in the fact that bisphosphonate therapy suppresses bone resorption, and thus results in a decrease in TGF‐β release, which then inhibits the proliferation of MM cells by suppressing IL‐6 secretion, because IL‐6 promotes the proliferation of MM cells 5, 8, 27, 28.

Osteoblast suppression occurring in MM patients facilitates the growth of MM cells and bone loss; therefore, bortezomib, which can enhance osteoblast proliferation, is suggested as a potential therapeutic intervention for MM. In its simulation here, bortezomib therapy is indeed shown to be effective in the management of MM‐induced bone disease through its action to decrease the viability of MM cells while limiting MM‐induced bone destruction. The inhibition of MM cells by the bortezomib therapy agrees with the experimental finding that increased osteoblast proliferation is able to reduce tumour burden in MM patients 23, 29, 30. Thus, the stimulation of osteoblast activity with therapies, such as bortezomib, can inhibit or even stop bone destruction as well as the tumour burden, and is an effective therapy for MM patients.

Because TGF‐β can indirectly promote the progression of MM cells, its inhibition is also suggested as a possible treatment for MM‐induced bone disease 5. However, although the model simulation indicates that this approach can lead to a decrease in MM cell viability, the MM‐induced bone loss is not inhibited and can become worse. In addition to the effect of TGF‐β on osteoblasts, TGF‐β can also inhibit osteoclasts by promoting their apoptosis 31, and TGF‐β inhibition would unavoidably lead to an increase in osteoclast activity and resultant bone resorption. This explains the increased bone loss resulting from TGF‐β inhibition. Therefore the MM management through inhibition of TGF‐β treatment does not appear to be effective using our modelling techniques.

The relationships between the treatment parameters (F.Bi, F.Bo and F.T) and equivalent drug dosages are currently not known; hence simulations to determine the sensitivity of the results to the values have been undertaken. The different treatment options can be clearly seen to work in different ways, and notwithstanding the uncertainty in the parameter values, all three simulations show qualitative agreement with the available clinical data, providing some degree of confidence in the model. Clearly further work and quantitative clinical data is required to confirm the parameter values and validate the model and its use in this application. Until then, it is hoped that this paper can serve as a virtual evaluation tool, which can be used to suggest new therapies or combinations of therapies and to explore the possible effectiveness of new therapeutic approaches before embarking on expensive clinical trials.