Osteoprotegerin in Bone Metastases: Mathematical Solution to the Puzzle

Bone is a common site for cancer metastasis. To create space for their growth, cancer cells stimulate bone resorbing osteoclasts. Cytokine RANKL is a key osteoclast activator, while osteoprotegerin (OPG) is a RANKL decoy receptor and an inhibitor of osteoclastogenesis. Consistently, systemic application of OPG decreases metastatic tumor burden in bone. However, OPG produced locally by cancer cells was shown to enhance osteolysis and tumor growth. We propose that OPG produced by cancer cells causes a local reduction in RANKL levels, inducing a steeper RANKL gradient away from the tumor and towards the bone tissue, resulting in faster resorption and tumor expansion. We tested this hypothesis using a mathematical model of nonlinear partial differential equations describing the spatial dynamics of OPG, RANKL, PTHrP, osteoclasts, tumor and bone mass. We demonstrate that at lower expression rates, tumor-derived OPG enhances the chemotactic RANKL gradient and osteolysis, whereas at higher expression rates OPG broadly inhibits RANKL and decreases osteolysis and tumor burden. Moreover, tumor expression of a soluble mediator inducing RANKL in the host tissue, such as PTHrP, is important for correct orientation of the RANKL gradient. A meta-analysis of OPG, RANKL and PTHrP expression in normal prostate, carcinoma and metastatic tissues demonstrated an increase in expression of OPG, but not RANKL, in metastatic prostate cancer, and positive correlation between OPG and PTHrP in metastatic prostate cancer. The proposed mechanism highlights the importance of the spatial distribution of receptors, decoys and ligands, and can be applied to other systems involving regulation of spatially anisotropic processes.


.1 Previous Work
The modeling part of this work is based on a previously developed model of bone homeostasis introduced in [1,2]. This model describes the spatio-temporal evolution of osteoclasts and osteoblasts as well as their interactions through the RANK/RANKL/OPG pathway, and can be summarized as follows Featured in (1) are the following state variables: osteoclasts (u 1 ), osteoblasts (u 2 ), the RANKL field (φ R ), the OPG field (φ O ) and the bone mass (z), which is ρ B in the current model. The fields y i are the active cell populations and differ from the total cell populations u i by an additive constant (see discussion in main text). Note that θ is the Heaviside function (θ(x) = 1 for x > 0 and θ(x) = 0 otherwise), and y 2,τ is a delay term defined by Among the various pathways involved in the remodelling process, only the RANK/RANKL/OPG pathway is modeled explicitly, the other pathways are captured by the nonlinearities in the u i equations. We discuss now the improvements and modifications of (1) that lead to the models used in Scenarios 1-5 in this study.

Osteoclast stimulation by RANKL
RANKL activates and stimulates osteoclasts by binding to RANK receptors on their surface. In (1), this ligand-receptor interaction is captured in the last term of the osteoclast equation as well as the third term in the RANKL equation, Note that we have two different reaction rates k 1 , k 2 to account for the reversibility of the RANK-RANKL binding: since a ligand can contribute to the stimulation of a cell without getting permanently bound to the receptor, we generally have k 2 < k 1 . The dependence on θ(y 1 )u 1 in (3) ensures that only osteoclasts in the vicinity of the remodeling front of active cells are stimulated by RANK-RANKL interactions. Even though this choice lead to satisfying results in previous studies, we chose to replace it as θ(y 1 )u 1 → y 1 . The reason for this choice is twofold: first, the dependence on θ(y 1 )u 1 implied that precursor cells responded to RANKL when close to the remodeling front, but not when further away from it. In absence of experimental evidence for such a functional distinction, the choice y 1 seems more meaningful. Second, from a mathematical point of view the use of the Heaviside function is expected to decrease the regularity of the solutions and hence by replacing it with a linear dependence as ∼ y 1 , the regularity of the solution is increased, and in particular the problem becomes more tractable from a numerical perspective.
Another issue with the RANK-RANKL term (3) is the proportionality constant k 1 : in fact, active osteoclasts are at all times attached to the bone surface and this implies that the activation of osteoclastic activity should vanish as the local bone density ρ B goes to zero. Despite these observations, the choice of keeping k 1 constant throughout this study (except for Scenario 3) can be justified a posteriori by observing that the RANKL gradients are consistently guiding the active osteoclasts away from previously resorbed spots (see e.g. Figure 8-A in main text). In Scenario 3 however, the RANKL gradient is reversed and guides osteoclasts towards the already resorbed centre of the domain (see Figure 6 in main text). Therefore, we modeled the linear dependence k 1 ∼ ρ B explicitly to avoid physiologically unfeasible situations as discussed in [2].

Chemotaxis
Chemotactic movement of osteoclasts has been documented in several studies, see e.g. [3], as well as [4] and [5]. Since distressed osteocytes have been shown to express RANKL [6,7], it seems plausible that there is a RANKL gradient away from damaged sites (e.g. microfractures). Interestingly though, there need not be a predefined RANKL gradient in the tissue: as illustrated by Figure 3-A in the main text, which starts with a homogenous RANKL field, the binding of RANKL to RANK receptors is sufficient to create a local gradient, which in turn drives the osteoclast movement. Furthermore, as pointed out in Scenario 4, there need not even be a sufficient RANKL background level at time t = 0: the production of PTHrP by the tumor is sufficient to induce the necessary RANKL for osteoclast stimulation and chemotaxis. In conclusion, the outcome of the simulations is largely independent of a pre-existing RANKL gradient in the tissue, and relies only on the well-documented properties of osteoclastic stimulation and attraction by RANKL.

Porous Diffusion
In the homeostasis model (1) we assumed φ R to obey porous diffusion dynamics, i.e. we allowed for exponents R > 1. The main difference between regular diffusion ( = 1) and porous diffusion ( > 1) becomes apparent when considering the model equation on R, The dynamics of (4) with = 1 have infinite propagation speed, i.e. for any t > 0 the solution φ will have infinite support, even if the initial field was compactly supported. On the other hand, solving (4) with > 1, compactly supported initial conditions will remain so for all t > 0, see [8]. Even though the choice of porous diffusion seems more meaningful for a mainly membrane bound cytokine such as RANKL, when we repeated all the simulations in this study for R = 1.5, 2, 2.5, we observed no qualitative differences to the case = 1. This observation, combined with the numerical advantage of working with a linear diffusion term, lead us to set R = 1 in this study.

Finite Osteoclast Speed
As observed in Figure 8A, tumour-derived OPG leads to an increase in osteoclast migration speed: the distance traveled by the remodeling front during 90 days monotonically increases with increasing OPG production τ O . The present model does not account for an upper limit as to how fast osteoclasts can move across the trabecula. Osteoclasts have been shown to travel at a speed of ∼ 10µm per day in physiological settings, and 5-10 times faster in cancer patients [9,10]. Since the maximal distance travelled by the cutting cone in our simulations is 7 mm in 90 days ( Figure 5C in main text), the numerical experiments reported in this study are within the range of realistic osteoclast speeds.

Undirected Osteoclast Motility
Many chemokines are known to induce random cell motility in addition to directional cell migration. To our knowledge, there is no experimental evidence that RANK-RANKL binding leads to an increase in undirected motility of osteoclasts. Nevertheless, we wanted to make sure that the nonlinear relationship Distance  Figure  5C of the main text are repeated, this time accounting for undirected motility of osteoclasts. The simulation is repeated for different values of initial RANKL φ 0 R , and different levels of OPG production by cancer cells, τ O . After 90 days, the following quantities are shown: distance traveled by osteoclasts (Distance), total number of active osteoclasts (OC), and total tumor mass (Tumor).
of tumor burden and OPG production in Figures 5C and 8B of the main text would also be observed if diffusive movement of osteoclasts did occur. In order to account for diffusive movement of osteoclasts, we replaced the osteoclast equation in system (6) of the main text by The new parameter σ u represents the diffusive motility of osteoclasts. In absence of experimental estimates, we used the tuning strategy for free parameters (as described in Section 3 below) to determine σ u . We set σ u = 10 −3 mm 2 d −1 since we found that higher values of σ u lead to osteoclast populations that moved too fast and did not stay compactly supported. We next assessed how osteoclast movement and tumor growth were affected by the changes in OPG production by tumor cells, and compared predictions of the model with osteoclast motility to the model without osteoclast motility (presented in Figure 5 in the main text). We found that accounting for osteoclast diffusion leads to a qualitatively similar behavior: the tumor burden remains a bell-shaped function of the OPG production rate τ O (Figure 1 of SI Text 1).
In view of these observations, we decided to avoid unnecessary complexity and uncertainty in our model, and set σ u = 0 for all simulations in the main text.

Reversibility of RANKL-OPG binding
We assumed that the binding of RANKL and OPG is irreversible, since to our knowledge, there is no experimental evidence that suggests otherwise. Nevertheless, to make sure that the possible existence of reversibility would not lead to qualitatively different results, we re-examined the simulations of Scenario 2, but this time accounting for reversibility of the [RANKL/OPG] binding. We introduced a new state variable φ RO and added the following reaction-diffusion equation to system (6) of the main text, where k 4 (units are d −1 ) is the dissociation rate, and σ RO is the diffusion constant for the complex. We assume that the complex remains membrane-bound, and therefore we set σ RO = σ R , the diffusion rate of membrane-bound RANKL (see Section 3). Furthermore, conservation of mass dictates the addition of the term +k 4 φ RO to the RANKL and OPG equations, respectively. Fixing an intermediate background RANKL level of φ 0 R = 2, we repeated Scenario 2 of the main text for three different values of k 4 = 0, k 4 = 0.1, and k 4 = 0.3, see Figure 2 of SI Text 1. The dissociation of the [RANKL/OPG] complex leads to a slight decrease of the distance traveled by the cutting cone, but it also results in higher osteoclast numbers in the cutting cone. Importantly, the qualitative impact of low to intermediate OPG production on tumor burden remains unaltered (compare to Figure 5C in main text).

Boundary Conditions
In [1,2], the simulations were performed on numerical domains that were large enough to avoid interactions between the bone remodeling unit and the domain boundaries. Therefore, Dirichlet boundary conditions in conjunction with spatial finite difference discretisations yielded good results. However, for the current implementation of the model we used periodic boundary conditions. This choice facilitated the use of spectral collocation methods for the discretisation of the Laplacian. In order to make sure that the size of the domain did not affect the simulations, we first solved the equations on the original domain, then doubled the domain size while keeping the mesh size constant, and eventually verified that the relative difference between the two solutions was negligible.

Fractional Step Method
The model equations of Scenarios 1-5 share a common feature: they all involve multiple time scales and hence suffer from stiffness, see e.g. [13]. In fact, the reaction and decay rates for the cytokines are of the order of an hour, whereas the chemotactic motion of the remodeling front has a time scale of about a week. If one were to perform a combined time-stepping for all state variables, the stiffness would strongly suggest the use of an implicit method. However, since most terms in the equations are nonlinear, an implicit time-stepping would be very expensive. Furthermore the chemotactic term, essentially hyperbolic in nature, is best tackled with an explicit method. These considerations lead us to solve the model equations with a fractional step method as proposed in [14]. Given an evolution equation with advective (A), diffusive (D) and reactive (R) contributions, the fractional step method consists of a sequence of iterative integrations. More precisely, given the value q n := q(t n ), the method integrates the solution over ∆t in three steps: 1. Advection. Solve ∂ t q = A(q) with initial datum q n over ∆t to obtain q * 2. Diffusion. Solve ∂ t q = D(q) with initial datum q * over ∆t to obtain q * * 3. Reaction. Solve ∂ t q = R(q) with initial datum q * * over ∆t to obtain q n+1 .
The main advantage of this approach is that for each contribution one can choose the best suited numerical method. For the advection step we used a second order centred difference scheme to discretize the spatial derivatives of the chemotactic term in the osteoclast equation. The time stepping between t and t + ∆t was performed by means of the Matlab built-in adaptive Runge-Kutta 45 solver ode45. Regarding the diffusion step we followed [14] and implemented TR-BDF2, an implicit Runge-Kutta method of second order. As for the spatial discretization of the Laplacian, we implemented a standard spectral collocation method, see e.g. [15]. Finally, the reaction step was integrated with ode45, similarly to the advection step. Even though the reaction part is mildly stiff for some Scenarios, ode45 performed well. To ensure that the employed numerical schemes were convergent, we performed convergence studies.

Parameter estimation
Similarly to the approach described in Section 4.1 of [2], we distinguish between fixed and free parameters. Fixed parameters are estimated based on experimental findings, whereas the remaining -a priori unknown -free parameters are determined according to the following tuning strategy. While keeping the fixed parameters unchanged, the free ones are tuned within a reasonable numerical range until the solution matches the following in vivo observations: the remodeling front of active osteoclasts stays spatially well-confined and moves across the trabecula at a roughly constant speed of 20 − 40 µm/day. In the remainder of this section, we implement this strategy and discuss individual parameters. Except for the parameters that vary from scenario to scenario, the values of fixed and free parameters are summarized in (7) and (8), respectively.
In a first step we consider the osteoclast (u), RANKL (φ R ) bone (ρ B ) and tumor (ρ T ) fields of Scenario 1, see equation (6) in the main text (ignoring the OPG and PTHrP fields). The parameters α, β and g, corresponding to the internal dynamics of the osteoclast population, have already been determined in [2] and are considered fixed. Since RANKL is assumed to be membrane bound, the experimentally unknown effective diffusion rate σ R should be at least one order of magnitude smaller than the diffusion rates of soluble OPG and PTHrP (whose values are discussed below). Performing simulations with σ R more than one order of magnitude smaller than σ O and σ P leads to numerical instabilities as is to be expected for an advection-dominated reaction-diffusion-advection system. In consequence, we choose σ R to be one order of magnitude smaller than σ O and σ P . If new experimental evidence would reveal that σ R is much smaller than assumed, our numerics would have to be reviewed. The remaining (free) parameters in Scenario 1, namely k 1 , k 2 , λ, ζ and k B could not be matched to any experimental data and hence we employed the tuning strategy outlined above.
For the additional fields appearing in Scenarios 2-5, namely OPG and PTHrP, as well as their interactions with the RANKL and bone fields, we have the following fixed parameters. First, the diffusion rates of OPG and PTHrP can be estimated from related experimental findings: in [16], effective diffusion rates for Verteporfin in subcutaneous and orthotopic tumours were measured to be 0.88 and 1.59 µm 2 s −1 respectively. Assuming that the diffusion rate σ of a molecule of mass M scales as σ ∼ M −1/3 , see [17], we can use these rates together with the molecular weights of OPG as a dimer (∼ 120 kDa [18]), PTHrP (∼ 18 kDa, [19]) and Verteporfin (∼ 0.7 kDa, www.drugbank.ca) to determine the range of possible effective diffusion rates of OPG and PTHrP in tumour tissue. For OPG we obtain a range of σ O ∼ 1.4 − 2.5 · 10 −2 mm 2 day −1 and for PTHrP σ P ∼ 2.6 − 4.7 · 10 −2 mm 2 day −1 . For the present study, we fix σ O = 1.6 · 10 −2 and σ P = 3 · 10 −2 . Note that the chosen scaling law σ ∼ M −1/3 applies to globular molecules, but may not be accurate for rod-like molecules. Since OPG and PTHrP have both globular and non-globular domains, a strict classification does not seem possible. Furthermore, we could not find a reliable source for the long and short axes of the molecules (which appear in the scaling relations for rodlike molecules), and hence we used the globular scaling relation to estimate the coefficients. Regarding the half-life of OPG, experimental data range between 10 minutes in rats [20] and 4 days in monkeys [21]. Converted into reaction rates this gives us a ballpark of k O ∈ [0.6, 110] day −1 . We decided to settle for an intermediate order of magnitude by choosing k O = 10 day −1 . Even though the half-life of PTHrP is expected to be of the same order of magnitude as the one of OPG, we could not find any corresponding experimental studies. Therefore, we treated k P together with the remaining parameters τ O , τ R , τ P and k 3 as free parameters and used the tuning strategy.
In summary, the fixed parameters are (we abbreviate day by d) α = 9.49 mm −1/2 d −1 β = 0.2 d −1 g = 0.5 σ R = 0.5 · 10 −2 mm 2 d −1 σ O = 1.6 · 10 −2 mm 2 d −1 σ P = 3 · 10 −2 mm 2 d −1 k O = 10 d −1 κ R = 1 d −1 , and the free parameters are A final note on parameter sensitivity is in order. Given the large number of parameters in the model, a systematic computational sensitivity analysis is not feasible. However, the sensitivity of the various parameters has been investigated in several steps. (A) Parameter sensitivity of the osteoclast-internal dynamics (production, apoptosis, autocrine stimulation) was assessed analytically in [22]. (B) Parameter sensitivity for the underlying spatio-temporal model of bone remodeling was studied by means of scaling arguments in [2]. (C) Sensitivity with respect to the parameters introduced in the current study (e.g. cytokine production rates by tumor, parameters in the PTHrP field) was assessed computationally over a wide range in parameter space.