Ten‐Year Simulation of the Effects of Denosumab on Bone Remodeling in Human Biopsies

ABSTRACT Postmenopausal osteoporosis is a disease manifesting in degradation of bone mass and microarchitecture, leading to weakening and increased risk of fracture. Clinical trials are an essential tool for evaluating new treatments and may provide further mechanistic understanding of their effects in vivo. However, the histomorphometry from clinical trials is limited to 2D images and reflects single time points. Biochemical markers of bone turnover give global insight into a drug's action, but not the local dynamics of the bone remodeling process and the cells involved. Additionally, comparative trials necessitate separate treatment groups, meaning only aggregated measures can be compared. In this study, in silico modeling based on histomorphometry and pharmacokinetic data was used to assess the effects of treatment versus control on μCT scans of the same biopsy samples over time, matching the changes in bone volume fraction observed in biopsies from denosumab and placebo groups through year 10 of the FREEDOM Extension trial. In the simulation, treatment decreased osteoclast number, which led to a modest increase in trabecular thickness and osteocyte stress shielding. Long‐term bone turnover suppression led to increased RANKL production, followed by a small increase in osteoclast number at the end of the 6‐month–dosing interval, especially at the end of the Extension study. Lack of treatment led to a significant loss of bone mass and structure. The study's results show how in silico models can generate predictions of denosumab cellular action over a 10‐year period, matching static and dynamic morphometric measures assessed in clinical biopsies. The use of in silico models with clinical trial data can be a method to gain further insight into fundamental bone biology and how treatments can perturb this. With rigorous validation, such models could be used for informing the design of clinical trials, such that the number of participants could be reduced to a minimum to show efficacy. © 2021 The Authors. JBMR Plus published by Wiley Periodicals LLC on behalf of American Society for Bone and Mineral Research.

In this study, a novel mechanoregulated fracture healing model is presented in which bone is modeled as a discrete entity, while the bone forming cells occupy this same spatial domain. The presented model combines paradigms of multiphysics simulation and agent-based modelling, creating a simulation in which discrete cells respond to their local mechanical and chemical environment by altering their substrate or producing more biomolecules.
A parameter study was performed using synthetic test cases and in vivo microcomputed tomography images of 20-week old female C57BL/6J mice, which had undergone an osteotomy. The effect of osteoblast polarization and the rate at which osteoid is produced was tested.
The model predicts that polarized osteoblasts are necessary for producing the porous hard callus observed during fracture healing, while the rate at which osteoid 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 98 is produced has a significant effect, not just in the rate of healing, but on the fate of the structure as well.

Introduction:
During the fracture healing process cells respond to a mixture of local and systemic signals which guide them in the regeneration of the damaged bone tissue (Marsell & Einhorn, 2011). Importantly, the signalling mechanisms involved travel across multiple scales to influence regeneration. At the organ scale the damaged bone is production, differentiation, mitosis, chemotaxis, chemokinesis, apoptosis and quiescence. Cells play a dual role by both sensing and influencing their environment. For example, cells influence their mechanical environments by modifying their chemical environment; osteoblasts and osteocytes in areas of low mechanical strain are known to produce less sclerostin (Robling et al., 2008) and increased amounts RANKL (Rubin et al., 2002). These circumstance lead to the inhibition of bone formation and the recruitment of osteoclasts which remove bone. Consequently this reduced structure will deform more under load and reduce the RANKL and sclerostin production. Conversely, a high local strain would decrease sclerostin production, increase OPG release and provoke osteoblasts to produce new bone tissue, which would stiffen the local structure reducing the local mechanical stimuli. Ultimately these two processes form a feedback loop which controls the microarchitecture of bone. The biochemicals produced by the cells also 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 99 propagate up to the organ level, with elevated levels of OPG and RANKL observed in both the fracture haematoma and the peripheral blood (Kӧttstorfer et al., 2014).
While pharmaceutical treatment through the systemic release of growth factors influence the amount of de-novo bone (Kruck et al., 2018;Gerstenfeld et al., 2009), and even the microarchitecture of the hard callus (Casanova et al., 2016). Changes to the microarchitecture are important as the structure determines the strength of the newly formed bone (Mehta et al., 2013). Thus, the microarchitecture of the callus, how it develops and how it is remodelled, is an important piece of the multiscale healing process. Incorporating this scale into simulations of fracture healing would improve the prognostic ability of fracture healing models.
A validated fracture healing model has the potential to reduce the need for animal testing when developing drugs and biomaterials. Current models can be categorised into three different groups: 1) mechanical models, in which the fracture healing process is governed entirely by mechanobiological rules. These models represent the tissues as a continuum. The majority of which have been developed for large animals (Ament & Hofer, 2000;Isaksson et al., 2008;Simon et al., 2011;Chen et al., 2009;Steiner et al., 2013;Vetter et al., 2011;Burke & Kelly, 2012). 2) Mechanics-free models, which are governed by rules derived from the cell biology of fracture healing. Current models translate the cell biology into a series of partial differential equations (PDEs), describing the movement, proliferation and differentiation of cells and their reactions with signal factors. As with the mechanical models the domains are modelled as continuums (Carlier et al., 2012;Geris et al., 2008;Bailon-Plaza & Van Der Meulen, 2001). 3) Combined models, which include a mixture mechanical regulation coupled with rules from the cell biology. These models have evolved from both the first two categories, combing elements of both. Geris et al. (2010) for example included mechanoregulatory elements from the work of Lacroix & Prendergast (2002) to their previous biological model. While mechanically regulated models have been expanded to include discrete cells which react to their local mechanical environment (Checa & Prendergast, 2009;Pérez & Prendergast, 2007). The use of discrete cells is becoming recognized retrospectively as a branch 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 100 of agent-based modelling (ABM) (Checa, 2018;Borgiani et al., 2017). Given the above, the current state of the field has two philosophies for modelling which are also compatible. However current models are limited in terms of resolution including only tissue level details as continuums and are either spatially low resolution or limited to two dimensions.
In this work, we build upon these existing models by combining a multiscale agentbased approach with mechanical and biochemical signalling. The discrete representation is extended to also include the tissue, which occupies the same spatial domain as the cells. Biochemicals are modelled as continua, but their receptors on cells are spatially discrete. Motivated by recent work showing the effect of WNT-signalling over the course of fracture healing (Kruck et al., 2018;Gerstenfeld et al., 2009;Naik et al., 2009;Kӧttstorfer et al., 2014), the simulation model includes the pathways related to OPG, RANKL, Sclerostin, VEGF and TGFβ.

Biochemical regulation of fracture healing
The fracture healing process is carried out by several actors. Cells from both the mesenchymal and hematopoietic families participate as shown in Mesenchymal stem cells are recruited from three sources: the bone marrow, surround tissue and the vasculature. It is currently believed that this process is driven by stromal cell-derived factor-1, BMPs and G-protein-coupled receptor CXCR-4 (Marsell & Einhorn, 2011). Once within the fracture callus, they are driven to differentiate into fibroblasts, osteoblasts and chondrocytes. These cells are responsible for building different types of extracellular matrix in order to stabilise the callus (Marsell & Einhorn, 2011). The hematopoietic stem cell (HSC) line is responsible for the formation of osteoclasts and immune cells. The osteoclasts are responsible for catabolic activity, removing unneeded material and remodelling the microstructure. Immune cells are vital actors in the inflammatory response, producing many cytokines which kick start the healing process, one of which is TGF-β which is responsible for upregulation of MSC proliferation (Stewart et al., 2010;Crane & Cao, 2014). There is also evidence that TGF-β down regulates the 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 101 proliferation of HSCs (Vaidya & Kale, 2015;Heino et al., 2002), and slows the differentiation of pre-osteoblasts into osteoblasts . The differentiation of MSCs into osteoblasts is regulated by the expression of BMP-2, which promotes differentiation into osteoblastic precursors, but inhibits the final differentiation into osteoblasts (Lee et al., 2003). There is growing evidence that the RANK-OPG axis is also important regulatory agents in fracture healing (Marsell & Einhorn, 2011;Kӧttstorfer et al., 2014). Gerstenfeld et al. (2009) have shown that the use of Denosumab, an anti-body for RANKL creates extremely dense callus tissue, in which mineralised cartilage is not resorbed, and the bone volume fraction (BV/TV) is significantly higher than the control group. As shown in Chapter 3.1, the resorption of mineralised tissue occurs in the third post-operative week, i.e. immediate after the formation of newly mineralised tissue. RANKL plays a crucial role in recruiting osteoclasts for this purpose. Another traditional bone remodelling protein is sclerostin. In a femoral midshaft osteotomy model, Kruck et al. (2018) showed that sclerostin antibodies enhanced bone formation, but could not prevent mechanically induced non-unions. Alzahrani et al. (2016) found that the callus which developed when treating mice with a sclerostin anti-body and a sclerostin knockout mouse 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 102 line was significantly larger than the control mice. As can be seen, the regulation of the fracture healing and bone remodelling processes have many commonalities.
Knowledge gained from the study of the mechanoregulation should be translatable from one to the other.

Modelling fracture healing
Agent-based models are starting to become widely adopted as a way of modelling biological systems. The implicit segregation of information within the model, as well as paradigms such as heterogeneous agents, make these models ideal for representing individual cells (Zhang et al., 2009). Elements of agent-based modeling have been used as parts of fracture healing models in the past, with some models incorporating cells as discrete individuals (Checa & Prendergast, 2009;Pérez & Prendergast, 2007;Kaul et al., 2015), and others including phenotype specific rules for tissue differentiation (Simon et al., 2011).
Solving chemical reactions: There are several methods for solving the reactions between ligands, receptors and anti-bodies at the continuum scale. A direct solution is to apply the law of mass action and generate an ordinary differential equation (Guldberg & Waage, 1879). This assumes the chemical species are well mixed and thus the rate of reaction is proportional to the species. The Hill-Langmuir equation, used in (Pivonka et al., 2008;Pastrama et al., 2018), is an equation which can be derived from the mass action laws. This describes the relative number of binding sites occupied, under the assumption of equilibrium between the forwards and backwards reactions. The Gillespie algorithm is a stochastic alternative to mass action. With the underlying assumption that the concentration of molecules is so low that the concentrations should be quantized as discrete molecules. In such situations the solution to mass actions represents the probability of a molecule being in a given location. On the other hand, the solution to the Gillespie algorithm is a single path through the probability space. For this reason, the Gillespie algorithm is most commonly used for reactions at sub-cellular resolutions (Gillespie, 1977). In this work we focused on solving the mass action equations directly. The focus on a direct solution rather than stochastic solution is 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 103 due to the scale at which the simulations are run. The Hill equation has several limitations regarding the choice of coefficients and the validity of the solution (Weiss, 1997). Coupling the non-linear reaction equation with the linear diffusion equation is another challenge. The method of lines (MOL) has been used by several authors in the field. MOL converts a partial differential equation (PDE) into a series of ordinary differential equations (ODE) which approximate the solution to the PDE. All equations are then solved iteratively with an ODE integrator. The issue which arises with this approach is when one part of the equation is stiff, the entire system must be solved with a small timestep. Such situations are common when solving coupled chemical reactions and diffusion (Sportisse, 2000). An alternative approach is operator splitting, where the reaction, diffusion and advection operations are applied sequentially (Strang, 1968). Operator splitting allows computationally intensive stiff operations to be solved independently of cheaper operations. It also allows specifically designed and optimized solvers to be used on each operation.
Development of mineralised tissue: During the healing process osteoid is deposited by osteoblasts, which is then mineralised to become either woven or dense bone. As this process comes to an end the osteoblasts have either become quiescent lining cells, undergone apoptosis, or have become osteocytes. The mechanism through which osteoblasts become embedded is not yet known. There are several possibilities, as described by Franz-Odendaal et al. (2006). In this study, we investigated two possibilities: 1) The osteoblasts are unpolarized, depositing matrix uniformly in their environment. 2) Osteoblasts are highly polarized depositing matrix directionally towards the surface upon which them sit.

Caveat: The combination of different modelling methods with a splash of biology
creates conflicts in the terminology used. For example, a "cell" for a biologist is clearly a biological cell, while for an engineer this could also be referring to a unit of volume on a computation grid. For this reason, the following nomenclature is defined: Cell, refers to a biological cell. Voxel, refers to a unit of volume, specifically cubic. Node, refers to a single computer in part of a HPC cluster. The developed simulation model combines paradigms from agent-based modelling, multiphysics and multi-scale modelling. In the following section, the ODD protocol is used to describe the model (Grimm et al., 2006). This splits the model description into three parts: 1) an overview which first gives context to the rest of the description by outlining the model purpose, introduces the elementary properties describing each cell, and describes the different length-scales and time-scales. 2) The design concepts, which introduce the expected behaviour of the cells and tissue as a collective. The cell specific rules are then introduced. 3) The details, this contains information important for implementation, but not necessary for understanding the model.

Purpose:
The purpose of this model is to provide an accurate mechanobiological representation of bone and in particular to provide an understanding of how complex micro-architectures can evolve and be maintained by individual cells responding to biochemical and mechanical signals.

State variables and scales:
The agents in the model represent individual cells in the bone tissue. Each agent occupies a voxel on a lattice uniquely and has a specific genotype. The genotype of the cell determines its behaviour, which ranges from movement within the lattice, adding or removing chemical signals or tissue, to differentiation into a different genotype. A complete list of cells in the model, as well as the binding sites and biochemicals the cells can produce, is listed in Table 4

.2: Schematic of simulation, A) An abstract piece of art showing healing callus with woven bone forming and vasculature infiltration. B) The coarse grid on which oxygen transport and consumption is solved. C) The lattice containing the cells, blood vessels and tissue. D) The lattices on which biomolecules diffuse. E) the internal states for an osteoblast.
is a direct linking between the lattices with the cell, with the osteoid and the mineral. Cells exist in the same space as these tissues and cannot co-inhabit a voxel 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 106 under normal circumstances. When a voxel has a greater than 50% occupation of osteoid the two following situations occur: 1) if the voxel is not occupied then no cell can enter it. 2) if the voxel is occupied by a cell this cell becomes embedded and can no longer move. The state variables for the voxels are the cell and all chemical concentrations, see Table 4.3. The lattice is enclosed with reflective boundary conditions, which represent the callus boundaries. Information can enter and exit the domain only via the vasculature cells, which act as a source/sink for most biomolecules and oxygen.

Process overview and scheduling
The simulation is designed to run over the course of 6 weeks, aligning with the measurements from the study of The structure of the bone is itself adapting throughout the simulation. While this is not modelled explicitly, the response of the cells to their chemo-mechanical environment directly leads to changes in the bone structure. The cells therefore use the local mechanical signal as a fitness function for their environment.

Interactions:
Interactions between cells occur through two possible ways: 1) Direct interactions can occur when cells are in contact (i.e. occupying adjacent voxels). 2) Indirect interactions occur when biochemicals are released from one cell and bind to another cell. An example of direct interaction are the osteoclast precursor cells which requires a certain number of osteoclastic neighbours before than can form or join an osteoclast cluster. An example of indirect interactions is osteocytes in regions of low strain recruit osteoclastic cells through releasing elevated levels of RANKL.

Stochasticity:
The migration of cells in the marrow is modelled as a stochastic process. Several surface cells will potentially migrate in random directions tangentially to the surface. The proliferation and apoptosis of cells is also a stochastic process. The vasculature development is also probabilistic. Finally, the differentiation of osteoblasts into pre-osteocytes can in some instances be related to stochastic behaviour of the cells.

Collectives:
Osteoclasts form clusters in which to function, while osteoblasts form layers. The formation of osteoclast clusters is partially stochastic, requiring at least 6 osteoclastic cells to be in contact. Their movement is stochastic but biased to move along the RANKL gradient. The formation of osteoblast layers is entirely stochastic, with osteoblast pre-cursors moving randomly along the surface until joining a layer, defined as a group of 3 or more osteoblasts.

.6 A mathematical framework for cell-substrate interaction and evolution
Several continuum level formulations exist for fracture healing and bone remodelling (Martin et al., 2017;Bailon-Plaza & Van Der Meulen, 2001). In this section, we describe a microscale mathematical framework for cell-substrate interactions. Firstly, the process occurs in a domain Ω = { ∈ 3 }. Within this domain we can define solid subdomains of bone and the osteoid as substrates upon which cells act. We suppose that these substrate domains Ω are a compact space bounded by ∂Ω , defined as: Where the substrate can be either osteoid or mineral. Within the domain Ω at time t there exists N(t) subdomains Ω , which are disconnected from all Ω such that Within the domain Ω there are cytokines which diffuse as follows: Where is the diffusivity of the cytokine in the serum. Each cell can act as a source for specific cytokines within the domain. The concentrations are updated using a per cell genotype update function: Where is a cell specific flux of cytokine across the cell membrane. When a new cytokine is produced by a cell and released omnidirectionally the change of the concentration in is: .
Where M is a cell specific function for producing the cytokine, and dS is the The directional production is thus: With this formulation a cell can directionally change the substrate, cells must do this while respecting the constraints of equation 2.
On the surface of the cells there can be binding sites which react with proteins outside of the cell domain in Ω. The nature of these reactions is assumed to be reversible with the following form Where the interaction distance is defined in Cells can proliferate and create new subdomains. This process is assumed to be random, with each cell having a specific probability based upon the internal states of the cell. The new sub domain is a copy of the existing one.

Lattice level
The primary simulation domain Ω is discretised as regular a 3-dimensional lattice.
The timesteps are discrete and the discrete voxels can be also be members of the subdomains for the cells or substrate. Changes to substrate concentrations on the lattice may redefine the substrate domain boundaries. As described in equation 2, voxels can either be occupied by a single cell or a tissue. There is an exception to this rule for osteocytes, which can be within the osteoid/mineral domain. Each site in the lattice can be identified with a coordinate vector The boundary conditions of the lattice are reflective such that for each axis a The occupied state can be either a freely motile cell , a vasculature cell , a surface attached cell , or a tissue , a summary of these sets can be found in Table   4.2. In the case that the voxel is occupied by tissue there are two state variables: The concentrations of both osteoid and mineral also define the substrate domains ( 1.

Cell motility
The freely motile cells within the marrow move according to a random walk in the Moore neighbourhood. The distance moved from the active voxel to the neighbour voxel is the random variable : For each Markov iteration a cell will move from ( ) to ( + ) if the location In the case an osteoblast is within the osteoid iso-surface it is considered trapped and differentiates into a pre-osteocyte.  HSCs move throughout the empty callus and marrow freely. In order for an HSC to differentiate into pre-osteoclasts, a sufficient level of RANK binding sites must be occupied, and the cell must have at least one voxel in the Moore neighbourhood with a mineralisation greater than 0.5. Pre-osteoclasts are not capable of resorbing bone. Instead pre-osteoclasts try to form clusters with other osteoclastic cells in order to become an osteoclast. They move along the iso-surface defined by > 0.5 . The movement is anisotropic, there is a 40% chance that they follow the gradient of RANKL, provided the gradient is non-zero, there is a 40% chance they move in a random direction along the surface, and finally a 20% chance they remain at their current location. As the surface is typically covered in cells, the pre-osteoclasts will swap position with cells in their desired location. Osteoclasts remain static until they are no longer attached to a surface voxel, once this occurs, they move as a pre-osteoclast. If an osteoclasts or pre-osteoclast cannot reattach to the surface it undergoes apoptosis.

Cell proliferation
In the model cells can duplicate themselves, simulating proliferation seen in nature.
This was a stochastic process in which random numbers were drawn for each cell.
If the number was above a given threshold the cell was copied to a neighbouring voxel. The chosen rates values are show in However, cell specific differentiation was used for it was assumed that in equilibrium the rates apoptosis and mitosis are equal, when MSC cells are stimulated by TGF-β the proliferation rate was increase, as were the HSCs with respect to binding of RANKL.

Cell mechanosensation
Each mechanically sensitive cell samples the local mechanical signal in the voxel it occupies. This in turn determines the amount of biochemicals and matrix the cells can produce. For both osteocytes and osteoblasts there is a linear relationship between the effective strain in the bone and the activation as described below: The activation is bounded between [0,1]. For osteoblasts there are two additional rule which takes precedent over this; if strains in the soft tissue are greater than = 0.3 then the activation is: And number of LRP5/6 binding sites occupied by sclerostin must be less than the number of unbound sites. Otherwise the cell become quiescent.
The activation is used as a scaling factor for the anabolic and anti-catabolic elements of the simulations, such as osteoid production, and OPG production. The complement of the activation is used for the catabolic and anti-anabolic elements such as sclerostin and RANKL. Osteoclasts and immune cells are a cell that once formed has a fixed activation of 1.0. A complete list of all cell types and their products can be found in Table 4.1.

Cell polarisation
The change in concentration of osteoid and mineral is divided into two phases. The first phase concerns the release of new chemicals for the cell and is as follows: The second phase determines the distribution of the newly produced material. This is described in equation 6. The following assumptions are made in order to solve these equations. Firstly, the size of the active surface is constant for all cell units of a given genotype. Secondly, that the size of the surface is dependent on the degree to which a cell is polarized rather than the local topology between the surface.
Finally, the infinitesimally small offset is taken to be the smallest unit of measurement on the lattice i.e. the voxel width. The orientation of the cell is determined as A normal distribution is creates at this location, which sums to unity in the Moore neighbourhood M of the cell. The change in concentration in a single voxel is then the weighted sum of all inputs from surrounding cells: Where τ is a voxel in the Moore neighbourhood, is a gaussian, and σ Is the standard distribution of the gaussian.

A multi-scale in silico model for fracture healing: Sensitivity of model parameters
The available cytokine was determined using equation 8. The discretization of which was as follows: Rather than a thin layer around the cell, the entire voxel the cell occupies is used as the region in which binding sites can react to the cytokine.
In addition to this simple case, the RANKL-RANK-OPG axis was also modelled, which can be described as a competitive reaction. With RANKL as the common target.

Mineralisation kinetics
While osteoblasts can produce osteoid, the mode through which it is mineralised is still not well understood. What is clear is that hydroxyapatite crystals condense at gaps within the collagen fibres reinforcing them (Nair et al., 2013). It is known that osteoblasts can produce vesicles containing hydroxyapatite nano-particles, which might acts as nuclei for the onset of mineralisation (Boonrungsiman et al., 2012).
However, such complex relationships were considered to be beyond the scope of this model. Instead, the mineralisation was considered a passive process. The amount of unbound osteoid represents the available sites for mineral to form, while the amount of available mineral in solute form is considered constant: 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters

122
The upper limit of the amount of mineral in a voxel is therefore the amount of osteoid within a voxel. In the proposed model a mineralisation rate of 0.125 per day was used. The mineral in solute was assumed to be a constant at 1.0.

Finite element simulations
The mechanical signal is calculated using the finite element method. For this purpose the solver ParOSol is used (Flaig, 2012). The boundary conditions and material properties are described in chapter 3.1. The mineral concentration was scaled linearly from the range [0,1] to a BMD in the range [0,800] mg HA/cm³. This BMD was converted to stiffness as described in chapter 2. For each voxel the effective strain was calculated. The conversion of effective strain to a mechanical signal was performed via gaussian dilation (Schulte et al., 2013). The soft tissue and hard tissue were treated independently.

Implementation of the lattice
The lattice was implemented as a C++ class which wrapped and provided access to the 3D data arrays. The grid containing the cells was an array of shared pointers to the parent class of the cells. This allowed cells of every genotype to be concurrently in the same array. AT the same time the use of shared pointers removed the need for active memory management and thus also simplified the code. All concentrations were stored as double precision floating points within the grid. The algebraic multigrids for the diffusion of molecules were calculated at model setup and stored in the parent class.

Implementation of numerical methods
The reaction-diffusion equations being solved are described in section 4.1.2.12 , and a brief overview of numerical methods used was discussed in section 0. In the presented model, the linear and non-linear operations are separated using second order Strang splitting (Strang, 1968). The diffusion was discretised with first order finite difference scheme with implicit Euler integration (BTCS). The resulting matrix was inverted using on a algebraic multi-grid solver from the open source AMGCL library (Demidov & Rossi, 2017).

Implementation of cell level parallelism
In the in vivo case every process is occurring simultaneously, in the in silico case this is not always possible due to issues with concurrent memory access known as race conditions. For example, two cells might try to enter the same voxel simultaneously, or to deposit osteoid in the same voxel. While trivial, these problems are dealt with two different ways. The movement of cells is discrete and sequential in time, in order to reduce simulation artefacts, the order in which cells are activated for movement is random. However, if this were purely random the race condition described above could occur if two processors moved two adjacent cells. To prevent this each processor receives two chunks of data to process (a first and second pass chunk), the chunks are spatially distributed such that data is only exclusively accessible for a single CPU during each pass. The processors can then 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 124 iterate randomly over the cells in each chunk with no risk of race conditions occurring. The same method is used for the production of osteoid, however random activation is not necessary, and thus is not used for performance reasons   An osteoblast polarization of = 0.7 was used for these simulations.

A multi-scale in silico model for fracture healing: Sensitivity of model parameters
127 The simulations ran on a Cray XC40/XC50 super computing system at the Swiss National Supercomputing Centre (CSCS, Lugano, Switzerland). The simulations used a single node and ran for approximately 2 hours.

In vivo micro-CT derived model
The main study used in vivo micro-CT images of the 0.85 mm group from chapter 3 as a basis. The model setup can be seen in  To investigate the effects of the discrete timestep for the cell behaviour an additional simulation was run on a single mouse. In this simulation, a larger timestep of 20 minutes for the cell movement, proliferation and matrix production was used.
For all in vivo simulations the osteoid production rate was 11,500 µm³/day, this was based upon results for the synthetic cubes. Osteoclasts had a fixed polarization of 0.7 and could resorb 11,500 µm³/day of both mineral and osteoid.
The simulations ran on a Cray XC40/XC50 super computing system at the Swiss National Supercomputing Centre (CSCS, Lugano, Switzerland). Each simulation used 6 nodes and ran for approximately 20 hours. The cubes were used to test the sensitivity of the rate of osteoid production in order to parametrise the larger in vivo simulations. All configurations resulted in "unions" in the third week, an example can be seen in Figure 4.8. The volume fraction of bone was proportional to the osteoid production rate as seen in   The structures deviate internally, for example the polarized osteoblasts create a structure which matches the native osteocyte density, with a overshoot at the time of fusion, Figure 4.16B. On the overhand, the unpolarized osteoblasts create a structure which at its peak has 60% more osteoblasts than the native structure, followed by a sharp decline at the onset of bridging in week 3, see Figure 4.16A. This is despite a higher concentration of RANKL in the fragment ROIs (Figure 4.18).
The RANKL concertation in the defect VOI rises for both parameter settings and promptly decreases post bridging. Week Highly Polarized High MSC seeding Unpolarized In vivo data     The remodelling is a crucial step in converting the structure from hard callus to cortical bone. The osteocyte density of the bone within defect VOIs (Figure 4.16) can explain the differences in the bone volume fraction. Unpolarized osteoblasts produce a larger surface of osteoid, leading to a higher recruitment of MSCs into osteoblasts. These results match the reports from Kaul et al. (2015), that only polarized osteoblasts can produce realistic distribution of osteocytes in tissue.
There is very little literature examining ultrastructural properties of bone, such as osteocyte density during the healing process. Casanova et al. (2016) found only a 4.1 A multi-scale in silico model for fracture healing: Sensitivity of model parameters 140 modest increase in the osteocyte density at three weeks post fracture when comparing the de-novo tissue to the existing cortical bone. Paralleling the trend observed for the simulations with polarized osteoblasts. It can be concluded that during fracture healing osteoblasts are polarized when depositing tissue, and that this is enough to embedded osteocytes at realistic densities.
Two initial conditions were tested, the marrow and periosteum were randomly seeded with either 1 or 2 million MSCs /mL. The high MSC density was shown to accelerate the early healing phase but result in an overall BV/TV not significantly different from the low MSC density group at the end of week 4. This can be explained by comparing Figure 4.10 and Figure 4.11, the initial amount of osteoid is higher for the highly polarized images, however there is an earlier fusion with locally denser tissue. Consequently, the entire defect is stress shielded which pauses the action of osteoblasts. This stress shield is evident in all groups, post bridging the sub-trochanteric region undergoes significant modelling because the load is transferred through the more mature tissue. There are two possible assumptions in the model which could be responsible for this unphysiological behaviour, 1) cells have no memory of their mechanical environment, a temporal averaging or leadlag filter would allow cells to slowly reduce their activity post bridging.
2) The rate of mineralisation could be reduced, this would delay the onset of mechanical bridging. However, the mineralisation rate already delays the mineralisation of an osteoid filled voxel by approximately an 8 day period. Other models have used time averaging of the mechanical stiffness matrix (Lacroix & Prendergast, 2002), while this was for numerical stability of their solver, it will have delayed the stress shielding post bridging. To summarise, further investigation is required regarding for how long cells remain stimulated post removal of the mechanical stimuli, and the rate at which osteoid mineralises in vivo.
The mathematical formulation defines the constraint regarding exclusive spatial occupation by either tissue or cells. This is the most significant change from existing models and is directly responsible for the formation of the microstructure. The mathematically formulation also describes the reaction between binding sites on there are other ways of solving this equation, firstly the simulations could be made more efficient by using a coarser grid for diffusion and via interpolation determine the number of molecules within the interaction distance. This would reduce the degrees of freedom but still satisfy equation 8. Secondly, the cell itself could be discretized into sub-cellular voxels using methods such as cellular Potts (Anderson & Rejniak, 2007). This would allow the cell to orient itself based upon binding site occupation and simulate phenomenon such as chemotaxis. Ultimately the presented the mathematical framework allows several modelling paradigms to be explored.
The sensitivity of the model to the timestep size is extreme. It is worth noting that halving the production of osteoid and maintaining the same timestep size in the synthetic cases had a far smaller effect compared to increasing the timestep between cell behaviours. The increase timestep size effectively halved the rate at which cells could move, this reduced the recruitment of MSCs to become osteoblasts, and the ability of pre-osteoclasts to cluster. Ultimately, this indicates that the fixed lattice for the cells is too fine for the 20 minute timesteps. Alternative methods such as particle based, or particle-in-cell based methods would allow subvoxel precision, continuous movement, and simultaneous movement of the cells.

Conclusion:
The