Bones have a number of different functions, but two of the most important ones are that bones provide structural support and physical protection of vital organs. Especially in elderly patients, bones can be fragile and frequently experience fractures. Although, bones have the fascinating and unique capacity to completely self-regenerate without leaving a scar; in some cases, such as in pathological fractures or those situations leading to large bone defects, bones fail to reach a complete healing. Delayed or incomplete bone regeneration represents a major clinical challenge. Nowadays, there is a need to understand the mechanisms involved in the bone healing process in order to develop treatment strategies which ensure the successful repair of the bone.

The unique process of bone healing is highly complex and dynamic and spans many different time and length scales. The process of healing can be divided into five overlapping phases: hematoma formation (pro-inflammation) phase, anti-inflammatory phase, soft callus formation (proliferative) phase, hard callus formation (maturing or modeling) phase, and remodeling phase (Schmidt-Bleek et al., 2014). During all these phases, many processes at the intracellular, cellular and tissue level scales are coordinated and interact to achieve bone restoration at the organ scale. At the intracellular level, a complex array of signaling molecules interacts and gives rise to the activation of specific genes that, ultimately, determine cell function. At the cellular level, cells proliferate, migrate, differentiate, and synthesize extracellular matrix. At the tissue level, bone, cartilage, and fibrous tissue are organized to provide the extracellular environment for the cells.

Over the last decades, experimental human and animal studies have provided detailed insights into the processes that occur during the bone healing response at individual length and time scales. For example, at the tissue level, it has been shown that bone healing is influenced by fixation stability (Schell et al., 2005), aging (Strube et al., 2008), fracture size (Claes et al., 1997), etc. At the cellular level, experimental studies have determined the contribution of the different cell phenotypes (immune cells, progenitor cells, etc.) to the healing response (Konnecke et al., 2014). At the intracellular level, the role of intracellular signaling pathways and their potential as therapeutic targets for bone regeneration have been investigated (Secreto et al., 2009). Despite great advances in the field at the individual length and time scales, the behavior of the process as a whole remains unclear.

The complexity involved in bone regeneration has brought computational modeling as a core discipline into regenerative research. In bone healing, such computer models aim to identify underlying rules driving bone regeneration cascades. First computer models of bone healing were focused on understanding how mechanical stability and/or loading influence the bone healing outcome. Using finite element techniques, those models aimed to determine the local mechanical signals (strains, pressure, fluid velocity, etc.) within the regenerating region and how they relate to the formation of different tissue phenotypes over the course of healing (Carter et al., 1988; Claes and Heigele, 1999). Although, they were a breakthrough in the field, they were mainly focused on the mechanical aspects of bone healing, ignoring or highly simplifying the biology of the process. As a result, they could not explain many of the experimental observations, such as non-healing in large bone defects or the effect of growth factor stimulation (Carlier et al., 2014; Ribeiro et al., 2015).

Over the years, computer models have progressively increased their complexity, by adding more and more biological details and giving rise to a new generation of computer models that try to understand mechanical and biological interactions occurring at the different time and length scales. One of the advantages of these more recent models lies in their ability not only to predict the bone healing outcome and how it might be influenced by, for example, fixation stability; but also to explain the mechanisms behind. Although mechanics can be used to predict the progression of bone healing in many situations, the mechanical conditions do not explain the bone healing process. Computer models that aim not only to predict but also to explain bone regeneration are based, at least to some extent, on the biological processes behind. Although some of these models have recently shown their potential to support the development of new treatment strategies for the promotion of bone repair (Carlier et al., 2014), most models have so far focused on understanding the mechanisms behind the bone healing response.

In what follows, we review how multiscale computer models have contributed to better understand the systems biology of the bone healing process. In the next section, we briefly describe the different modeling approaches that have been used to simulate bone regeneration. Thereafter, we group the models in three main categories based on the main biological processes under investigation: cellular activity (including migration, proliferation, differentiation, extracellular matrix production), growth factor production and effect, and angiogenesis.

Modeling approaches in the field of bone regeneration have mostly focused on finite element, partial differential equations (PDEs) and agent-based techniques (Table 1). Finite element analyses have been used to solve the mechanical equations, i.e., to determine the mechanical signals within the regeneration region. Many models have then coupled these mechanical signals to changes in biological parameters; for example the differentiation or migration of the cells (Lacroix et al., 2002; Gomez-Benito et al., 2005; Isaksson et al., 2008; Nagel and Kelly, 2010; Checa et al., 2011). Following an iterative approach, these models simulate how changes in tissue patterning influence mechanical signals within the healing region and how these signals further influence tissue formation (Lacroix et al., 2002; Isaksson et al., 2006; Checa et al., 2011). Systems of PDEs have been used to simulate temporal and spatial changes in cell and tissue density or protein concentrations over the course of healing (Bailon-Plaza and van der Meulen, 2003; Gomez-Benito et al., 2005; Geris et al., 2010b). This modeling approach focuses on the processes occurring at the tissue and cell population (i.e., cell density) level. Models that take into account the processes at the cellular level are based on agent-based approaches. These models employ experimentally derived computational-coded rules to define the behavior of individual cells (Byrne et al., 2011; Checa et al., 2011). These models simulate the bone regeneration process at the cellular scale and aim to understand how cellular behavior gives rise to tissue formation pattern and bone regeneration.

Tissue patterning over the course of bone healing is an emergent event that results from multiple processes taking place at the lower length and time scales. For example, tissue patterning is determined by the number of cells present in the healing region, their location, migration, and proliferation capacity, as well as by the amount of extracellular matrix produced by the different cell phenotypes. Bailon-Plaza et al. were the first to explicitly introduce biological factors into a computer model of the bone healing process, including simulation of cell migration, proliferation and differentiation, as well as production and resorption of corresponding tissues (Bailon-Plaza and van der Meulen, 2001). Using PDEs, they determined spatial and temporal changes in the density of mesenchymal stem cells, osteoblasts and chondrocytes, as well as changes in bone and cartilage tissue densities inside the callus. Although, the model presented some major limitations, such as a highly simplified callus geometry, they were able to predict early periosteal bone formation and endochondral ossification in the external callus according to experimental observations. Thereafter, Bailon-Plaza et al. (Bailon-Plaza and van der Meulen, 2003) and others (Gomez-Benito et al., 2005; Isaksson et al., 2006, 2008; Garcia-Aznar et al., 2007; Geris et al., 2008, 2010a,b; Burke and Kelly, 2012) adapted and extended this model to investigate different aspects of the bone healing process, such as the effect of mechanical loading on the bone healing process. Bailon-Plaza et al. (Bailon-Plaza and van der Meulen, 2003) simulated the stimulatory or inhibitory effect of mechanical signals on the ossification process according to Claes et al. (Claes and Heigele, 1999). With this model, they were able to explain the beneficial and adverse effects of moderate and excessive loading, respectively, as well as the negative effects of delaying mechanical stimulation of rigidly fixed calluses; as observed experimentally. Geris et al. investigated the occurrence of non-unions due to mechanical over-loading (Geris et al., 2010b). In their model, they showed that the mechanical regulation of both angiogenesis and osteogenesis was able to predict overload-induced non-union, confirming the hypotheses of experimental studies investigating the interconnection between angiogenesis and osteogenesis. All these models based on PDEs make use of diffusion equations to simulate cellular invasion within the callus. Several models have used this approach to investigate the effect of the origin of the cells on the healing pattern (Lacroix et al., 2002; Vetter et al., 2012). Comparing model predictions with histological data, Vettel et al. showed that the source of progenitor cells needs to be stronger at the periosteum than at the marrow space in order to achieve early periosteal and delayed endosteal bone formation, as seen in bone healing experiments in sheep (Vetter et al., 2012).

Several computer models have used agent-based approaches to simulate the behavior of individual cells and to investigate their contribution to tissue patterning over the course of healing (Byrne et al., 2011; Checa et al., 2011; Carlier et al., 2012; Borgiani et al., 2015; O'Reilly et al., 2016). Using this approach, Checa et al. (2011) and Borgiani et al. (2015) showed that differences in the healing pattern (e.g., endosteal, periosteal, time to bringing) of different species (sheep, rat, and mice) can be explained by differences in the regulation of cellular activity by mechanical signals during the course of healing.

Statistical methods have been used to investigate the robustness of the models and to identify the influence of model parameters on model outcome. Using a design of experiments approach, Isaksson et al. determined the rate of bone tissue formation and cartilage degradation as key players for the prediction of uneventful bone fracture healing (Isaksson et al., 2008). Using the same approach, Carlier et al. identified the impaired endochondral ossification process and increased infiltration of fibroblastic cells as key contributors to the degree of severity of congenital pseudarthrosis of the tibia (Carlier et al., 2016).

Growth factors play a key role in the bone healing process by modulating cellular activity and serving as communication between the cells. A single growth factor may have effects on multiple cell types and induce different functions, resulting in a cascade of highly complex interactions. Several computer models have aimed to understand the role of growth factors on cellular activity, and thereof on tissue formation over the course of healing (Bailon-Plaza and van der Meulen, 2001; Moore et al., 2014; Ribeiro et al., 2015). Bailon-Plaza et al. assumed two initial sources of an osteogenic and a chondrogenic growth factor located at the periosteum and the hematoma, respectively; while keeping the duration of the release and the magnitude as parameters. They predicted that the duration of the release in the hematoma needs to be sufficiently long to initiate the healing response. They found an agreement with experimental data in the time point of peak osteogenic growth factor concentration (Bailon-Plaza and van der Meulen, 2001).

The chemotactic potential of vascular endothelial growth factors on the angiogenic process during bone healing has also been investigated using a computer model (Geris et al., 2008). This model suggested that endothelial cells are attracted toward gradients of vascular endothelial growth factor secreted by chondrocytes and osteoblast and that these gradients are ensured by the consumption of the growth factor by the endothelial cells. Based on these assumptions, the model was able to simulate compromised healing conditions due to impaired angiogenesis. However, the model did not take into account the influence of the mechanical environment on growth factor production or its effect on vascular growth.

Moore et al. explicitly included the effect of mechanics on the production of bone morphogenetic protein 2 (BMP-2) by progenitor cells located in the periosteum of a critical sized femoral defect in sheep stabilized with an intramedullary nail (Moore et al., 2014). Assuming that BMP-2 in turn regulates cell proliferation and differentiation, they were able to predict experimentally observed tissue formation over the course of healing.

The potential of growth factor BMP-2 incorporation to promote cellular activity and enhance tissue regeneration in the context of large bone defects has been shown not only in a clinic setting (Schmidt-Bleek et al., 2016), but also using computer modeling techniques (Ribeiro et al., 2015). Ribeiro et al. collected quantitative experimental data from the literature on the effect of BMP-2 on cellular activity (proliferation, migration, differentiation, extracellular matrix production, etc.) and implemented them into a model to simulate bone healing in a large defect stimulated with BMP-2 (Ribeiro et al., 2015). The authors showed a good agreement with experimental data on the amount of bone tissue formed over the course of healing; however they were not able to predict endosteal bone formation.

Angiogenesis, defined by the formation of new capillaries, plays a key role in bone repair. Several computational studies have been developed to investigate the relative role of revascularization on bone healing progression (Geris et al., 2008; Chen et al., 2009; Peiffer et al., 2011; Simon et al., 2011; Burke and Kelly, 2012; O'Reilly et al., 2016). These models have seen a clear evolution toward a multiscale approach with the aim to better reflect and understand the biology of the process. While early models simulated vascular growth as a simple diffusion process (Geris et al., 2008; Chen et al., 2009; Burke and Kelly, 2012), more recent models use a discrete approach in order to simulate the complex structure of the newly formed capillary network (Carlier et al., 2015; O'Reilly et al., 2016).

Diffusion based models have investigated the influence of angiogenic growth factor production or callus size on the establishment of the vascular network and in turn on successful bone healing (Geris et al., 2008; Chen et al., 2009). Computer models at the cellular level have focused on understanding how progenitor cell fate might be affected by the altered conditions present in fractures with disrupted vasculature (O'Reilly et al., 2016) or in genetically modified mouse animal models (Burke et al., 2013). The MOSAIC model, developed by Carlier and her group, is the only multiscale computer model of bone healing that incorporates processes at the intracellular scale (Carlier et al., 2012). It implements an intracellular module that includes Dll4/Notch1 signaling to determine tip cell selection. The model was able to simulate the salt and pepper pattern seen for cell fates, an increased tip cell density due to the loss of Dll4 and an excessive number of tip cells in high VEGF environments (Carlier et al., 2012).

Computer models can provide great insights into the complexity of the bone regeneration process, especially on the interactions between the different length and time scales. Despite significant advances in the field of computer modeling of bone healing and its clear trend toward a multiscale approach to provide a systems biology overview of the process, many limitations remain. In what follows, we describe some of the remaining challenges and give insights into possibilities for future work.

Despite advances in modeling being rather rapid, the present computer models are still unable to capture many of the mechanical and biological interactions that occur across and within the different time and length scales. More importantly, in most cases, they lack prediction power in compromised repair conditions. Old age (Bak and Andreassen, 1989), immune compromise (Claes et al., 2012) or genetic disorders (El Khassawna et al., 2012) have been shown to negatively affect the bone healing response; however little is known about the mechanisms behind these phenomena. From a computational point of view, very little has been done to understand how these factors influence the bone healing cascade. In silico studies show a great potential to contribute to the understanding of bone healing in compromised conditions and the consequences of alterations in cellular function on the bone healing outcome.

EB and SC drafted the manuscript; all authors read and revised the manuscript and approved its content.

This study was supported by the German Research Foundation (Deutsche Forschungsgemeinschaft: DU298/14-1; CH 1123/4-1).