Towards in silico Models of the Inflammatory Response in Bone Fracture Healing

In silico modeling is a powerful strategy to investigate the biological events occurring at tissue, cellular and subcellular level during bone fracture healing. However, most current models do not consider the impact of the inflammatory response on the later stages of bone repair. Indeed, as initiator of the healing process, this early phase can alter the regenerative outcome: if the inflammatory response is too strongly down- or upregulated, the fracture can result in a non-union. This review covers the fundamental information on fracture healing, in silico modeling and experimental validation. It starts with a description of the biology of fracture healing, paying particular attention to the inflammatory phase and its cellular and subcellular components. We then discuss the current state-of-the-art regarding in silico models of the immune response in different tissues as well as the bone regeneration process at the later stages of fracture healing. Combining the aforementioned biological and computational state-of-the-art, continuous, discrete and hybrid modeling technologies are discussed in light of their suitability to capture adequately the multiscale course of the inflammatory phase and its overall role in the healing outcome. Both in the establishment of models as in their validation step, experimental data is required. Hence, this review provides an overview of the different in vitro and in vivo set-ups that can be used to quantify cell- and tissue-scale properties and provide necessary input for model credibility assessment. In conclusion, this review aims to provide hands-on guidance for scientists interested in building in silico models as an additional tool to investigate the critical role of the inflammatory phase in bone regeneration.


INTRODUCTION
Bone healing is a complex, well-coordinated process that starts autonomously when a bone fracture occurs. Bone fractures are one of the most common injuries and their incidence in Europe is expected to increase by 23% over the coming decade due to ageing, as average life expectancy rises (Borgström et al., 2020). Owing to the bone tissue characteristics, successful healing is usually achieved within weeks (Marsh, 1998). However, up to 10% of bone fractures result in delayed healing or non-union (Zura et al., 2016). This risk rate is influenced by anatomical location, fracture severity and host factors such as age, smoking or the presence of comorbidities (Zura et al., 2016;Mills et al., 2017;Stewart, 2019). Current treatment options to prevent or cure these incidences present many drawbacks. Autologous bone grafting remains the gold standard procedure to treat nonunions, but this technique has limitations such as significant donor site morbidity and limited volume of available tissue (Pape et al., 2010). Alternative approaches to support the healing process, such as bone tissue engineering strategies, are still being tested in clinical trials or under development (Amini et al., 2012;Lammens et al., 2012;Papantoniou et al., 2021). These approaches mainly target the skeletal system and the repair phase of fracture healing, whereas recent findings have demonstrated that the skeletal and immune system are closely interacting through a carefully coordinated cross-talk between inflammatory and bone forming cells. Hence, inflammatory cells, such as macrophages, are believed to play a critical, but yet incompletely understood role in bone healing (Schlundt et al., 2018;Pajarinen et al., 2019).
In the last two decades, computational modeling has developed into a powerful technique to complement and reinforce traditional in vitro and in vivo experimentation, as it can provide an integrated view of the many events happening during the bone healing process and hence lead to a deeper understanding of said process. Moreover, computational models aim to reduce animal experimentation, although in vivo studies are often still required for validation purposes. Despite the importance of the inflammatory phase in bone fracture healing (Könnecke et al., 2014;Loi et al., 2016;Schmidt-Bleek et al., 2016), most computational models of bone regeneration focus on the repair phase, ignoring inflammation and its impact on the regenerative outcome. Therefore, modeling inflammation is a necessary inclusion in the current state-of-the-art, as it will allow to elucidate the mechanisms regulating the early phase of bone healing and their effect on the final regenerative outcome. This inclusion is not only needed, but it also starts now to be possible, as an increasing amount of experimental data regarding the inflammatory response is becoming available in the literature (Könnecke et al., 2014;Kovach et al., 2015;Schlundt et al., 2018;Wagar et al., 2018;Maruyama et al., 2020).
To date, despite the advances in the experimental work on the inflammatory phase of bone healing, no in silico models exist that capture the spatiotemporal dynamics of the process. In this review, we bring together the necessary components required to build a validated computational model able to predict the inflammatory response in bone healing and study the interaction of this phase with the subsequent phases of the healing process. First, Section 2 describes the overall bone fracture healing process with a strong focus on the inflammatory phase. Then, Section 3 revisits the current computational models describing the immune and the skeletal system responses after an injury. Next, Section 4 presents an overview of the experimental techniques that can be used throughout the development of computational models, from calibration to validation. Finally, the concurrence of biology, computational methods and experimental validation is discussed in Section 5. Taken together, this review aims to provide the necessary information and tools to build in silico models, which can provide an additional perspective to study the critical role of the inflammatory phase in bone regeneration.

THE BIOLOGY OF BONE FRACTURE HEALING
Bones support the body, enable its mobility and protect vital organs. Moreover, bones produce hematopoietic cells and contribute to mineral storage within the bone marrow. Bone tissue is highly dynamic: bones adapt themselves to changes in the body, accommodating mechanical and biological requirements, and are constantly renewed in a process of remodeling. However, when stress and compression forces overcome bone tissue tolerance, bone fracture occurs (Oryan et al., 2015) and the process of fracture healing starts.
Bone fracture healing is described in Subsection 2.1. Subsection 2.2 focuses in more detail on the inflammatory response during fracture healing, paying special attention to cellular activity, cytokines and mechano-regulation.

Bone Fracture Healing Process
Bone can regenerate autonomously without fibrous scar formation after most cases of injury or fracture, eventually restoring its original state. This healing capacity is orchestrated by the complex fracture healing process, which involves multiple different cell types and is regulated by several biochemical, physical and mechanical factors (Einhorn, 1998). Depending on the mechanical stability of the fracture, direct or indirect healing will occur. Direct or primary fracture healing leads to restoration of the bone through a remodeling process. However, primary fracture healing is rather exceptional as it requires complete stability at the fracture site (Marsell and Einhorn, 2011), which is typically not achieved (Perren, 2002;Harwood et al., 2010;Claes et al., 2012). On the contrary, indirect or secondary fracture healing, the most common form of fracture healing (Marsell and Einhorn, 2011), is stimulated by interfragmentary motion (Harwood et al., 2010;Claes et al., 2012). In secondary fracture healing, bone repair advances via a multi-staged process involving both intramembranous and endochondral ossification (Loi et al., 2016), in which bone is formed directly from mesenchymal tissue or from intermediate cartilaginous tissue, respectively. However, high interfragmentary motion inhibits bone healing progression (Claes et al., 2012), resulting in compromised healing.
The classic phases of secondary fracture healing are inflammation, repair and remodeling ( Figure 1). This simple classification is further elaborated in the contemporary literature, where additional overlapping substages have been proposed: hematoma formation, acute inflammation, granulation tissue formation, angiogenesis, fibrous tissue formation, fibrocartilage, soft callus development, cartilage mineralization, hard callus development, and, finally, remodeling (Kolar et al., 2010;Loi et al., 2016). Following Figure 1, these key events are briefly described below.
Immediately after the trauma, the fracture hematoma is formed due to the blood vessels disruption, which triggers the blood coagulation cascade, thus creating a fibrin network. This fibrin network serves as provisional extracellular matrix for the influx of inflammatory cells as well as the progenitor cells from the periosteum and the bone marrow (Kolar et al., 2010;Loi et al., 2016).
Although this phase of bone healing is mostly defined as the invasion of inflammatory cells, the hematoma also contains immune cells present in the blood released from the disrupted vessels (Kolar et al., 2010). The healing process is then initiated with the activation of neutrophils, monocytes and macrophages (Kolar et al., 2010), leading to an acute inflammation reaction and the release of growth factors and cytokines. The initial fracture hematoma (Grundnes and Reikerås, 1993;Kolar et al., 2010) and subsequent inflammatory response (Bastian et al., 2011;Claes et al., 2012;Hoff et al., 2016) are critical for fracture healing (Schlundt et al., 2015). The hematoma is cleared in several days by the action of macrophages, which remove the fibrin matrix and necrotic cells at the bone ends via phagocytosis (Loi et al., 2016). A hypoxic environment remains within the fracture site, since neovasculature has not been developed yet.
At the end of the inflammatory phase, granulation tissue replaces the hematoma fibrin network due to the recruitment and proliferation of skeletal progenitor cells (SPCs) and fibroblasts (Harwood et al., 2010;Marsell and Einhorn, 2011). Granulation tissue favors angiogenesis, which is the formation of new blood vessels from pre-existing ones (Carano and Filvaroff, 2003). The vascularization process of the fracture site is promoted with interfragmentary motion during the early stages of fracture healing (Claes et al., 2012) and enhanced with angiogenic factors, such as fibroblast growth factor (FGF), platelet-derived growth factor (PDGF) or vascular endothelial growth factor (VEGF) (Barnes et al., 1999;Carano and Filvaroff, 2003;Tsiridis et al., 2007). Meanwhile, the hypoxic environment in the central region of the fracture site induces the differentiation of SPCs into chondrocytes (Tsiridis et al., 2007;Claes et al., 2012), starting the repair phase. Chondrocytes produce cartilage to connect the fractured bone ends, forming a soft callus that wraps the fracture gap. The soft callus provides initial mechanical stability and serves as scaffold for endochondral ossification during the repair phase (Harwood et al., 2010;Marsell and Einhorn, 2011;Loi et al., 2016). At the same time, SPCs differentiate into osteoblasts in the periosteal region away from the fracture site, hence creating woven bone via intramembranous ossification (Malizos and Papatheodorou, 2005;Claes et al., 2012;Loi et al., 2016). Both ossification events are regulated by growth factors, such as bone morphogenetic protein (BMP) and transforming growth factor beta (TGF-β), which control proliferation, differentiation and apoptosis of both chondrocytes and osteoblasts (Barnes et al., 1999;Tsiridis et al., 2007). As soft callus chondrocytes proliferate, they become hypertrophic and secrete VEGF, generating the rightful environment to attract blood vessels. Hypertrophic chondrocytes will finally undergo apoptosis and blood vessels will recruit progenitor cells that will differentiate into osteoblasts, leading to cartilage mineralization and generating the hard callus (Marsell and Einhorn, 2011). The formation of the hard callus entails the end of the repair phase of bone healing, leaving a solid and mechanically rigid fracture site, which has been revascularized and repopulated with bone cells. This stage is reached within several weeks or even months after the trauma (Loi et al., 2016) and generates the mechanobiological conditions to initiate the process of bone remodeling.
The remodeling phase is the final stage of the bone healing process and takes years to complete. Bone remodeling involves the resorption of immature woven bone and underlying cartilage matrix by osteoclasts, replacing these tissues with lamellar bone, as well as the decay of osteoblasts, undergoing apoptosis, or their maturation and embedding into the bone matrix as osteocytes. The cellular functions of osteoclasts and osteoblasts are regulated by cytokines such as receptor activator of nuclear factor kappa B ligand (RANKL) and osteoprotegerin (OPG). While RANKL promotes cell activation, differentiation and survival, OPG inhibits cell activation and induces apoptosis (Steeve et al., 2004;Tsiridis et al., 2007); leading to bone remodeling being controlled by the RANKL/OPG ratio (Steeve et al., 2004). The remodeling process establishes the osteon structure and Haversian system of the bone, restoring the bone's original shape, strength and stability (Oryan et al., 2015;Loi et al., 2016).

Inflammatory Response to Bone Fracture
The inflammatory response is the immediate reaction to a trauma that starts when pathogenic agents enter the body due to a wound, generating undesired living conditions for the injured organism (Ward and Lentsch, 1999;Bastian et al., 2011;Loi et al., 2016). The inflammatory response consists of an automated cascade of signals that activates the innate immune system to contrast the invasion (Osuka et al., 2014). Cells of the innate immune system can directly attack the pathogen or trigger a second wave of signaling by releasing specific factors that can support the response (Stoecklein et al., 2012). Environmental conditions, such as swelling and temperature increase, are generated within the injured zone due to external attack encouraging a quicker inflammatory response (Evans et al., 2015;Loi et al., 2016).
The inflammatory response has a primary role on the overall bone healing process, and immune restricted patients are more prone to experience impaired healing (Hoff et al., 2017). For example, cytokines and growth factors released from macrophages recruit SPCs, promoting their capacity to colonize the fracture zone and differentiate, thus progressing the healing (Loi et al., 2016). Slower completion of bone fracture healing is observed in case of reduced influx of macrophages (Alexander et al., 2011;Schlundt et al., 2018). However, a continuously activated inflammatory response may incur a chronic state, which is also detrimental to successful healing (Osta et al., 2014). This chronic inflammatory fate was observed in numerous cases of delayed bone healing, where the prolonged exposition of the healing tissue to cytotoxic T cells extended the pro-inflammatory stage to the detriment of a fast and successful healing (Schmidt-Bleek et al., 2012). Adequate treatment of bone fracture healing should therefore generate a balanced response from the inflammatory stage. While it is known that the antiinflammatory environment generates the conditions for a successful repair phase (Godwin et al., 2017), the prolonged use of nonsteroidal anti-inflammatory drugs was observed experimentally to alter the healing process (Lisowska et al., 2018). Inflammation involves a large number of agents that cooperate at different time and length scales to guarantee an adequate response. In the following subsections, we will describe the characteristics and functions of the principal immune cells and cytokines that are involved in this process.

Cells of the Immune System
The cells involved in the inflammatory response can be divided into two groups according to their belonging to the innate or adaptive immune system (Medzhitov and Janeway, 1997). The cells of the innate immune system, which include monocytes, macrophages, neutrophils, natural killer cells and dendritic cells, constantly monitor the organism and provide the first response to the pathogens (Medzhitov and Janeway, 2000;Bouchery and Harris, 2019). The adaptive immune system guarantees the second pathogen-specific reaction and is mainly regulated by the migration of T and B lymphocytes, also referred to as T and B cells, within the infected region. This response is not immediate and requires more time to process and enter into action. However, the adaptive immune system can keep a copy of the antigen to accelerate the response in case of a future attack from the same pathogen. A full characterization of immune cells and their role in the inflammatory response is beyond the scope of the review and it is already well described elsewhere (Mosser and Edwards, 2008;Chaplin, 2010).
Cells of both the innate and adaptive immune system are present in the fracture site during the inflammatory stage of bone healing (Baht et al., 2018). For example, circulating neutrophils and monocytes migrate to the healing region in the first hours after the injury (Hoff et al., 2016;Kovtun et al., 2016). Neutrophils are the first cells to be recruited in the healing region to promote the formation of the fibrin thrombus to stabilize the fracture (Bastian et al., 2016). Monocytes circulate within the bloodstream, ready to extravasate from the capillaries to the surrounding tissues when the inflammatory response is triggered (Maslin et al., 2005). In the bone healing scenario, monocytes are also recruited from the bone marrow and they invade the fracture region to clean it from debris and to upregulate the pro-inflammatory response (Soltan et al., 2012). Once in the fracture site, monocytes will turn into adherent cells and differentiate into macrophages. The macrophages present within the fracture gap in the early stage of the inflammatory response are activated by the pro-inflammatory environment. Traditionally, macrophages were described to be activated into two states, named classically (M1) and alternatively (M2) activated, depending on whether they promote or inhibit the inflammatory response. Macrophage activation within the two states is fundamental for the right course of the inflammatory phase of bone healing. M1 macrophages regulate the initial proinflammatory response and clean up the region from dead cells and pathogenic agents through phagocytosis (Mescher, 2017). Additionally, they promote the recruitment of other proinflammatory cells through the secretion of specific cytokines (Gu et al., 2017). At the end of the inflammatory phase, the macrophages differentiate into M2 macrophages, which downregulate the inflammatory response to create the right environment for the following repair phase. The release of adequate anti-inflammatory cytokines provokes the recruitment of repair cells, such as SPCs and fibroblasts, which will start the rebuild of the fractured bone. The importance of the role of macrophages as initiators of the repair phase has been shown by Schlundt et al. (2018), who observed altered endochondral ossification in cases where the macrophages within the fracture site were depleted. Recently, experimental work has shown that macrophage activation is rather like a spectrum than a two-state system, with specific signatures depending on the location within the spectrum (Mosser and Edwards, 2008;Harasymowicz et al., 2021).
The adaptive response in fracture healing starts when T cells sense the molecular signals released from the cells of the innate immune system within the injury region. If the injury is characterized by infection from an external pathogen, specific antibodies are produced and released from B cells to accelerate the neutralization of the threat. The adaptive response in fracture healing follows a two-wave dynamics as it is observed to peak both after the fracture and later during the cartilage revascularization (Könnecke et al., 2014). Contrasting ideas are reported on the role of the adaptive immune response on the overall bone regeneration process: while fracture healing is observed to be accelerated when the adaptive reaction is suppressed (Toben et al., 2011), the positive effect of certain categories of T cells on bone regeneration has been also reported (Sadtler et al., 2016;Jahn and Weidinger, 2017). Furthermore, T and B cells promote the differentiation and recruitment of osteoclasts, which shift the balance of the later remodeling phase to favor bone resorption over formation (Manabe et al., 2001;Gillespie, 2007).

Cytokines
Besides cells, the inflammatory response is mediated at subcellular level by molecular signals called cytokines. Cytokines are small proteins released from the inflammatory cells to regulate the inflammatory response, thus playing a role in the correct development of the early stages of bone healing. According to their ability to enhance or inhibit the inflammatory response, cytokines are typically divided into pro-inflammatory or anti-inflammatory, respectively.
Pro-inflammatory cytokines, such as interleukin 1 and 6 (IL-1, IL-6) and tumor necrosis factor alpha (TNF-α), are observed to peak their expression in the healing region within the first 24 h post-injury (Dimitriou et al., 2005). At subcellular level, TNF-α regulates the correct development of the inflammatory phase and, furthermore, it enhances the recruitment of SPCs to initiate the subsequent repair stage (Karnes et al., 2015).
Anti-inflammatory cytokines, which include different interleukins such as IL-4, IL-10, IL-11 and IL-13, are released to reduce inflammation when the first wave of pro-inflammatory response is over. Anti-inflammatory cytokines downregulate the inflammatory response and prevent chronic inflammation, which would be detrimental to fracture healing (Zhang and Yao, 2019).

Mechano-Regulation
Although cells and cytokines are the major biological regulators of the inflammatory response in bone fracture healing, the micromovement in the interfragmentary region also regulates bone healing at cellular level. The early stage of bone healing is particularly sensitive to changes in mechano-stimulation, hence establishing an adequate mechanical environment at the injury site is necessary from the beginning of the inflammatory phase (Klein et al., 2003).
Monocytes, for example, express a stronger pro-inflammatory response under shear or compressive loading (Fahy et al., 2019). Mechano-regulation also affects the behavior of macrophages during the inflammatory phase as tissue stiffness influences their activation status, shape, mobility and phagocytic capacity (McWhorter et al., 2013;Adams et al., 2019;Jain et al., 2019;Gruber and Leifer, 2020). Elongation of macrophages under influence of mechanical loading induces anti-inflammatory activation and initiates the repair phase in the healing process (McWhorter et al., 2013). Thus, adequate fracture mechanical support is decisive to shape the macrophages and move from the inflammatory to the repair phase (Ballotta et al., 2014).

IN SILICO MODELING
With the term in silico, scientists refer to the wide field of research that benefits from the use of computer modeling and simulation to investigate intricate and complex systems. This approach is becoming established in the biomedical field, as an additional resource to obtain a detailed understanding of the organism or its individual components. The flexibility provided by the computational approach favors the unveiling of aspects and insights that would be otherwise challenging to monitor experimentally. For this reason, the use of in silico clinical trials in all stages of the research and development pipeline has progressively gained more attention in the last decades (Viceconti et al., 2016). One of the more recent applications of in silico models is the execution of in silico clinical trials. In this context, the use of in silico models (e.g. through the use of Monte Carlo methods or the Bayesian approach) allows to quantify the parametric uncertainty in large data sets obtained from the results of the computational simulations. This is one way in which the effect of population variation can be captured in silico. To date, this approach is used by researchers to investigate the mechanisms behind neuron activation (Marder and Taylor, 2011), action potential stimuli in cardiac cells (Britton et al., 2013;Lawson et al., 2018), or mechanical properties of skeletal muscles (Sierra et al., 2015), among others.
There is a wide range of in silico models available in the literature to investigate different aspects of bone regeneration in silico (see Doblaré et al., 2004;Giorgi et al., 2016;Borgiani et al., 2017;Ghiasi et al., 2017 for recent reviews), but most of them focus only on the repair and remodeling phases, thus ignoring the inflammation phase. At the same time, many approaches modeling immune and inflammatory responses in other tissues have been presented in the last two decades (see Section 3.2). The inflammatory response is already a complex process to simulate and, certainly, additional complexity arises when it is included in the bone healing model, as interactions and processes at different biological levels (tissue, cellular and subcellular) have to be considered. Besides, there are different computational approaches (continuous, discrete or hybrid) to model these different levels over their relevant time and length scales. The main characteristics of each approach are discussed in Subsection 3.1 and illustrated in Figure 2A, in order to elucidate which approach might be best to use depending on the biological goal of the research. Next, an overview is provided of the most relevant models to investigate the inflammatory response (Subsection 3.2) and the bone regeneration process (Subsection 3.3). In addition, Subsection 3.4 introduces the only in silico model that is, to the best of the authors' knowledge, currently investigating the inflammatory response in bone fracture healing.

In silico Approaches to Model Biological Processes
There are different mathematical approaches to model a given biological situation: deterministic or stochastic, continuous or discrete over time and length scales, phenomenological or mechanistic. The appropriate approach is determined by the question that needs to be answered, the context of use, the available data and the computational resources. Even then, various models can produce qualitatively similar behavior (Anderson and Chaplain, 1998;Alber et al., 2006). Since the interest of an in silico model of the inflammatory response in bone fracture healing lies in understanding the biological events happening within the fracture region, we will focus here on models that incorporate physiological processes. To build such type of model, a common framework in mathematical biology is the so-called compartmental model.
A compartmental model is a system with different compartments and transitions between them ( Figure 2B). In Figure 2B, specific biological entities (chemical factors, cells and extracellular matrices) have been assigned to one compartment depending on their type and can interact with entities in other compartments by transitions equipped with rates. As a result of interacting compartments, a coupled system of conservation equations is derived, in which each compartment is represented by one equation. Transitions between compartments represent biological processes such as migration, differentiation or apoptosis. Rates often follow the law of mass action and are modeled using rate formulations such as Michaelis-Menten kinetics or the Hill function. If the biological FIGURE 2 | In silico approaches to model the bone healing process and the inflammatory response. (A) Overview of in silico techniques to describe biological processes and predict their different outcomes. The choice of the in silico model depends on the research goal. Continuous models are often used to describe general dynamics at tissue and cellular scales, such as bone mechanics, in which different tissue matrices interplay (figure adapted from Wang and Yang, 2018). Discrete models are mostly used to represent individual behavior at (sub)cellular scales, such as the immune response, which comprises a high number of cells and cytokines. The hybrid approach combines the advantages of both continuous and discrete techniques, providing comprehensive multiscale models that allow to investigate, for instance, sprouting angiogenesis during the bone regeneration process (figure obtained with the model described in Carlier et al., 2016).  1 | Overview of the biological agents and processes in bone fracture healing and the way they can be captured in silico in continuous-time or discrete-time models. In keeping with the description in Section 3.1, the spatial scale (if present) is mentioned. Examples are provided of experimental setups to use during the calibration phase of in silico models. Units of quantitative parameters that can be extracted from experiments are in squared brackets. For those experiments that lead to qualitative observations, this is mentioned explicitly. entities migrate either randomly or directed up to a gradient (such as chemotaxis or haptotaxis), a diffusion term is considered in the equation, providing a solution as function of time and space. Additionally, density-dependent models include growth dynamics by using e.g. a logistic-growth function (Murray, 1989). More details about common principles to model biological processes (in bone fracture healing) can be found in Table 1.
Compartmental models can be translated into deterministic or stochastic models, using continuous-or discrete-time approaches. In most cases, the solution of the resulting system of equations is impossible to solve and represent using analytical techniques and hence it is approximated with numerical methods. Continuous-time models use differential equations to describe mechano-biological processes. Differential equation models, whether ordinary (ODE), delay (DDE), partial (PDE) or stochastic (SDE), imply a continuous overlap of generations (Murray, 1989), thus describing the chronological time of biological phenomena. ODEs often describe how spatialaverage biological entities change over time, simulating e.g. inflammatory responses (Kumar et al., 2004;Reynolds et al., 2006;Vodovotz et al., 2006;Trejo et al., 2019) or the bone healing process (Trejo et al., 2019;Lo et al., 2020) at tissue and cellular levels, and individual intracellular dynamics (Warrender et al., 2006;Peiffer et al., 2011) at subcellular level. DDEs model biological processes that not only depend on the current time, but also on an earlier time; representing e.g. hematopoiesis regulation (Mackey and Glass, 1977;Faria and Oliveira, 2020) or inflammatory responses (Nagaraja et al., 2014). PDEs describe the spatiotemporal evolution of biological entities, using e.g. reaction-diffusion equations to model bone healing within the fracture area (Bailón-Plaza and Van Der Meulen, 2001;Lacroix and Prendergast, 2002;Gómez-Benito et al., 2005;Isaksson et al., 2006;Geris et al., 2008 among others, see section 3.3) or angiogenesis (Olsen et al., 1997;Anderson and Chaplain, 1998). SDEs introduce random parameters in the model and are used to investigate e.g. bone remodeling (Sun et al., 2012;Jerez et al., 2018). One main advantage of continuous-time models is that they have been studied exhaustively in the last centuries, leading to many well-known numerical methods to determine their solutions, such as the finite element (FE) method (Zienkiewicz et al., 1977, used e.g. in Zysset et al., 2013Coquim et al., 2018), the finite difference method (Liszka andOrkisz, 1980, used e.g. in Nagatani et al., 2006) and the method of lines (Hundsdorfer andVerwer, 2003, used e.g. in Geris et al., 2008). Some main disadvantages are that they usually fail to capture heterogeneous behaviors (Van Dyke Parunak et al., 1998;Wilkinson, 2009) and that the incorporation of new biological aspects is often not trivial (Yates et al., 2001;Fachada et al., 2007).
Discrete-time models use difference equations to study smallscale biological processes at (sub)cellular levels. Difference equations do not consider overlap between successive generations as they are solved for each time increment, involving an inherent delay to register changes (Murray, 1989). Discrete models are characterized by a stochastic nature, allowing the introduction of probabilistic rules, such as Monte Carlo methods (Lux, 2018), to describe each biological entity with its own properties and not as part of a population (Alber et al., 2006). The most common discrete approaches in biomedicine are agentbased (AB) and cellular automata (CA) models. AB models simulate the behavior of agents that can evolve generation after generation by changing their spatial position and internal properties. AB models are typically used to investigate cellular dynamics in response to environmental conditions, finding many applications in 1 | (Continued) Overview of the biological agents and processes in bone fracture healing and the way they can be captured in silico in continuous-time or discrete-time models. In keeping with the description in Section 3.1, the spatial scale (if present) is mentioned. Examples are provided of experimental setups to use during the calibration phase of in silico models. Units of quantitative parameters that can be extracted from experiments are in squared brackets. For those experiments that lead to qualitative observations, this is mentioned explicitly.  (1998) Frontiers in Bioengineering and Biotechnology | www.frontiersin.org September 2021 | Volume 9 | Article 703725 immunology (Martínez et al., 2012;Shi et al., 2016;Pappalardo et al., 2020 among others, see section 3.2) and bone healing (Checa et al., 2011;OReilly et al., 2016;, 2019 among others, see section 3.3). CA models are a subgroup of AB models, hence they are often referred to as AB. However, while AB models are more focused on the single agent behavior to explore its impact on the overall scenario, the CA method is based on nearest-neighbor interactions governed by phenomenological rules (Anderson and Chaplain, 1998), meaning that interactions in CA models regard only neighbor regions. CA models are typically used to describe angiogenesis (Bentley et al., 2008;Peiffer et al., 2011;Carlier et al., 2012). One advantage of discrete-time modeling is its capacity to model each single element as an individual entity, allowing heterogeneous behaviors (Van Dyke Parunak et al., 1998;Wilkinson, 2009). Some main disadvantages are that the number of unknown parameters is usually high, and these rely crucially on the biological parameters obtained from experimental data (Murray, 1989). This often entails a reduction of precision and accuracy, resulting in model simplifications and approximations (Shi et al., 2016). Continuous and discrete models can complement each other in hybrid models. In hybrid models, continuous-and discretetime approaches are coupled through input/output variables to provide multiscale models describing mechano-biological processes at tissue, cellular and/or subcellular levels. The most common hybrid formulation in biomedicine couples a PDE system of reaction-diffusion equations with AB or CA models (Stéphanou and Volpert, 2016). The PDE system is often used to capture biomechanical stimuli and tissue mechanical properties (Checa et al., 2011;Zahedmanesh and Lally, 2012;Virgilio et al., 2015;Ceresa et al., 2018), which regulate the spatial distribution of cells modeled with an AB model. CA models are often used to describe angiogenesis, and coupled to PDE systems describing the spatiotemporal evolution of cells, tissue matrices and chemical factors (Peiffer et al., 2011;Carlier et al., 2012). It is also common to use continuous formulations to regulate the subcellular behavior of individual cells within a discrete model. For instance, ODEs regulating the intracellular behavior of endothelial cells (Peiffer et al., 2011;Carlier et al., 2012) or PDEs determining the molecular environment of individual cells (Warrender et al., 2006). The reader is referred to Stéphanou and Volpert (2016) for a review of hybrid modeling in biology.

Modeling the Inflammatory Response
Computational models describing the immune response, also known as computational immunology, can be broadly classified into two groups: those describing generic inflammatory responses after infection or trauma, and those simulating immune responses in specific tissues. Most approaches of the former group are continuous, whereas the latter are often simulated using discrete or hybrid models. Starting with the models describing generic inflammatory responses, Kumar et al. (2004) proposed a three-equation ODE model to simulate a simplified acute inflammatory response, able to predict healthy and negative outcomes. This model describes the relationships between the pathogen, which instigates the innate immune response, and early and late pro-inflammatory mediators (Kumar et al., 2004). Reynolds et al. (2006) considered three subsystems for different biological situations (non-specific local immune response, resting phagocytes and activated phagocytes) and merged them into a four-equation ODE model describing the generic acute inflammatory response to a pathogen. A bifurcation analysis of the model identified when the outcome was compromised depending on the administration of antiinflammatory mediators (Reynolds et al., 2006). In the same year, Vodovotz et al. (2006) introduced a more elaborate mathematical model to simulate a non-specific acute inflammatory response after trauma, infection or hemorrhagic shock. This ODE system described the dynamics of cells and cytokines and included the effect of tissue dysfunction, coagulation elements and blood pressure. In addition, it was the first model validated with animal and human experimental data . Almost a decade later, Nagaraja et al. (2014) presented a comprehensive mathematical model to represent the local inflammation process in a wound and characterize the indicators triggering chronic inflammation. This model was validated with experimental data and consisted of fifteen ODEs and one DDE: the ODEs described inflammatory cells, cytokines and growth factors, whereas the DDE represented monocyte differentiation into proinflammatory macrophages, driven by chemotaxis with a 12 h delay (Nagaraja et al., 2014).
Within in silico models of tissue-specific immune responses, the popularity of AB methods is clear (Fachada et al., 2007;Shi et al., 2016). Several AB models of the immune system can be found in the literature, together with a large variety of simulators to develop them. Some models are implemented in custom AB simulators, such as ImmSim Seiden and Celada, 1992) and UISS (Pappalardo et al., 2010), whereas other models use generic open-source simulators that allow for the implementation of AB models such as NetLogo (Mi et al., 2007;Brown et al., 2011;Pennisi et al., 2013;Shi et al., 2016). ImmSim was the first AB model and framework to simulate the immune system, focusing on the processing of antigens and their effects on the different cell types Seiden and Celada, 1992;Fachada et al., 2007). Mi et al. (2007) presented an AB model focused on the interrelation between inflammation and skin wound healing in a physical domain. Skin injury and the subsequent inflammatory response were simulated to examine the general healing progression in terms of cells and cytokines dynamics. Brown et al. (2011) described a model of inflammation simulating the response of macrophages and fibroblasts to particulate exposure in the lung, as well as their interactions within the simulated environment, such that cytokines production, tissue damage and collagen deposition are represented. Martínez et al. (2012) developed a model of macrophage action on endocrine pancreas, focused on modeling the activation of the innate immune system upon stimulation by necrotic or apoptotic cell death in the first step of type 1 diabetes autoimmune response. Wendelsdorf et al. (2012) designed the ENISI simulator to represent mucosal inflammatory and regulatory immune pathways in the gut. Shi et al. (2016) proposed an integrated-mathematical-AB model to simulate the hepatic inflammatory response to Salmonella infection in mouse, which might cause a severe immune response and result in sepsis. Pennisi et al. (2013) and Pappalardo et al. (2020) investigated the cause of chronic inflammation in relapsing remitting multiple sclerosis using AB techniques. This framework was further developed into the Universal Immune System Simulator (UISS), which is now also used to investigate immunotherapy in cancer (Gianì et al., 2018) and the development of vaccines for Tuberculosis (Russo et al., 2020b) and Sars-Cov2 (Russo et al., 2020a).

Modeling the Repair and Remodeling Phases in Bone Healing
The use of computer models to simulate bone healing can be dated back to Carter et al. (1988), who investigated in silico the role of intermittent stress on the revascularization and tissue differentiation processes in the initial stages of bone healing. In the following years, many other studies exploited the computational power to study the mechano-regulation of the bone healing process (Prendergast et al., 1997;Carter and Beaupré, 1998;Claes and Heigele, 1999;Bailón-Plaza and Van Der Meulen, 2003;Isaksson et al., 2006;Geris et al., 2010;Checa et al., 2011;Burke and Kelly, 2012;Vetter et al., 2012;Borgiani et al., 2015;Wang and Yang, 2018). Another common application of computational methods is the simulation of the revascularization process on bone healing to highlight the role of angiogenesis and relative oxygen supply on disrupted tissues (Geris et al., 2008;Peiffer et al., 2011;Carlier et al., 2012;Carlier et al., 2015;OReilly et al., 2016). Moreover, different in silico models have been developed to investigate critical healing therapeutic strategies, such as the use of bone graft with a scaffold support (Perier-Metz et al., 2020), the transplant of stem cells (Geris et al., 2010;Carlier et al., 2016) or the provision of exogenous growth factors (Moore et al., 2014;Ribeiro et al., 2015) within the healing region. Different biomechanical studies employed in silico approaches to evaluate the impact of fracture stabilization (Gómez-Benito et al., 2006), gap size (Gómez-Benito et al., 2005) and nature of mechanical stimuli (Epari et al., 2006;García-Aznar et al., 2007;Steiner et al., 2013).
Most of the aforementioned studies use FE analyses to reproduce the mechanical environment (e.g. stress/strain distribution, tissue mechanical properties, bone density) within the injury. However, to date, many studies in this field started to additionally employ AB models to acquire a different point of view on the investigation of the mechano-biological relationships driving bone fracture healing. The supporting AB models are commonly employed to simulate the dynamics of repair cells (Byrne et al., 2011;Checa et al., 2011;Borgiani et al., 2019Borgiani et al., , 2021 and angiogenesis (Peiffer et al., 2011;Carlier et al., 2012;OReilly et al., 2016). Carlier et al. (2012) developed the hybrid MOSAIC model to simulate sprouting angiogenesis in a discrete environment. The behavior of the discrete endothelial cells was regulated by their protein levels and their relationship with cells, tissue and growth factors present in the global continuous environment. The multiscale model from Checa et al. (2011) investigated the inter-species differences in bone fracture healing between small and large animals within a mechano-regulated environment. They used an AB model to simulate how the spatial distribution of specialized bone repair cells (microenvironment) is regulated according to the mechanical stimulus predicted with FE (macroenvironment). Multiscale in silico modeling is a successful approach to explore the bone healing process at the levels of tissues, cells and subcellular agents by simulating their response to mechanobiological stimuli.

First Model of the Inflammatory Response in Bone Healing
The in silico studies of the repair and remodeling phases reported in the previous section do not include the simulation of the early stages of bone fracture healing, thus ignoring the role of the inflammatory response. Inflammation is characterized by numerous actors whose role in the overall scenario is worthy to be investigated. However, due to its complex nature, the inflammatory response to bone injury has been rarely simulated with a computer model. To the best of the authors' knowledge, only one in silico model describing the inflammatory response in bone fracture healing has been reported in the literature. The model was first introduced by Kojouharov et al. (2017) and further updated within the same research group by Trejo et al. (2019). Kojouharov et al. (2017) developed an eight-equation ODE model to simulate the temporal dynamics of debris, cells, cytokines and tissues from the first hours post-fracture, capturing the interaction between biological elements acting at multiple levels. Debris removal was modeled with a constant rate depending on the debris and macrophage densities, while the macrophages density depended on migration and emigration rates. The concentration of pro-and anti-inflammatory cytokines was simulated using Hill functions to capture a saturation effect, which depended on the concentration of debris and macrophages and of SPCs, respectively. Finally, the dynamics of SPCs, osteoblasts, fibrocartilage and woven bone was described as in Bailón-Plaza and Van Der Meulen (2001). The model simulated the biological time-dynamics in different case scenarios, highlighting the influence of a controlled cytokine concentration level as treatment to obtain an overall successful Frontiers in Bioengineering and Biotechnology | www.frontiersin.org September 2021 | Volume 9 | Article 703725 healing. Moreover, the model was employed to propose cytokinebased treatment in challenging healing conditions. For example, the model showed faster acceleration when an optimized dose of anti-inflammatory cytokines was administered at the beginning of the healing process. Two years later, Trejo et al. (2019) incorporated two additional equations to simulate the distinction between classical and alternatively activated macrophages, and the ODE system was adapted accordingly ( Figure 2B). Other biological processes were upgraded as well, the most relevant being debris removal, modeled now by a Hill function to represent the saturation of phagocytosis by macrophages, and macrophages migration, described now by a logistic growth function. The updated model allowed to analyze the role of macrophage activation status in the inflammatory phase to generate a successful signaling cascade initiating the subsequent repair phase. The model endorsed macrophages as promoters of tissue production during healing, giving further merit of this enhancement to the alternatively activated ones (M2). However, no spatial distribution of the different biological agents was modeled, as only temporal evolutions were reported as results.
The spatiotemporal evolution is of upmost importance to further explore the dynamics of all the involved actors during the progress of the inflammatory response, as it represents the heterogeneous distributions within the region of interest. For example, one hypothesis would be that macrophage migration to the fracture zone will initially have a bigger impact in the peripheral area and less effect in the central area, generating different spatial dynamics in the healing process. Therefore, we believe that the next generation of in silico fracture healing models should include both temporal and spatial evolution of the densities and concentrations of the different biological agents related to the inflammation phase. Moreover, many experimental studies investigate the immune response nowadays, as presented in Section 4. The incorporation of the spatial description in silico would allow a stronger validation of the future computational models investigating bone fracture healing from the initial inflammatory response to the later remodeling phase.

EXPERIMENTAL VALIDATION OF IN SILICO MODELS
Experimental techniques are continuously evolving to study the inflammatory response on multiple scales, ranging from microscale in vitro systems to large in vivo animal models. These results provide important information also in view of validating the predictive capacity of bone healing in silico models. Each modeling technique has its unique advantages and provides essential information about the inflammatory process ( Figure 3A). In vitro models allow the culture of human cells in a controlled environment outside of living organisms, although they are poorly suited for long-term studies. Moreover, in vitro Frontiers in Bioengineering and Biotechnology | www.frontiersin.org September 2021 | Volume 9 | Article 703725 models may fail in recapitulating a clinically relevant environment due to the absence of all factors present in vivo (Boussommier-Calleja et al., 2016), which motivates the use of animal models. The resemblance of the human biological environment is the reason why in vivo models are an absolute requirement for translational studies of human immunology (Wagar et al., 2018). However, biological mechanisms may differ between animal models and humans (Mestas and Hughes, 2004). Hereafter, both traditional and advanced in vitro and in vivo systems to model the inflammatory response in bone healing are discussed in Subsection 4.1 and Subsection 4.2, respectively. In Subsection 4.3 we describe different assays to extract both qualitative and quantitative data for the validation of in silico predictive models.

Source of Inflammatory Cells
Human blood is the most frequently used source of immune cells for in vitro experiments since peripheral blood samples are easy to obtain. Immune cells with a single nucleus can be isolated from the whole peripheral blood by density centrifugation (Dagur and McCoy, 2015). These cells, named peripheral blood mononuclear cells (PBMCs), are a heterogeneous cell population mainly composed of lymphocytes and monocytes. Lymphoid cells account for 85% of all human PBMCs and consist of T cells (∼60%), B cells (∼10%) and natural killer (NK) cells (∼10%). Monocytes constitute around 15% of the total PBMCs count, while other cell types, such as dendritic cells, are less than 1% (Bittersohl and Steimer, 2016).
In general, in vitro experiments study specific cellular functions and require the isolation of single cell types. While monocytes are traditionally isolated from the rest of PBMCs and differentiated into macrophages by cell adhesion to tissue culture plastic (Rios et al., 2017), every cell type in PBMCs can be separated by labeling with magnetic beads. Immunomagnetic cell separation consists of binding magnetic beads to cell surface antigens using specific antibodies (Plouffe et al., 2015). The characteristic surface molecules, named cluster of differentiation (CD) molecules, of each PBMC type are known: CD3 + for T cells, CD22 + for B cells, CD56 + /CD16 + for NK cells and CD14 + for monocytes (Bittersohl and Steimer, 2016).
Although human PBMCs are routinely used to answer fundamental questions about the immune cell functions, their choice has some drawbacks. First, immunomagnetic cell separation requires expensive reagents, such as antibodies. Secondly, in vitro results obtained from primary human cells are affected by natural immune variations between individuals, which is related to genetic variations, environmental exposure and aging (Patin et al., 2018). A variable immune response is crucial in the context of patient-specific models of the immune response, but it might obscure the effects related to the mechanisms under investigation. Therefore, in search of a stable phenotype, human cell lines are an attractive solution for many in vitro experiments aiming to validate in silico models. As an example, THP-1 is an established monocytic cell line isolated from a patient affected by acute monocytic leukaemia (Tsuchiya et al., 1980). Compared to human primary monocytes, THP-1 cells can be cultured in vitro for an almost indefinite time, while maintaining monocytic characteristics. On top of that, there is limited genetic variation between THP-1 cells, thus their phenotype is stable during culture. Nevertheless, the polarization profile of THP-1 cells does not coincide with the one of primary monocytes isolated from PBMCs. It is suggested to use THP-1 cells to validate in silico models involving phagocytosis and M1 activation (Shiratori et al., 2017).

Mono-Culture vs. Co-Culture
Single immune cell types have been extensively investigated in traditional mono-culture systems such as tissue culture plastics. Macrophages, for example, are routinely derived from monocytes and activated using standard activators, such as IL-4, IL-10, TGFβ, interferon gamma (IFN-γ) and lipopolysaccharide (LPS). Each standard activator, or their combination, is associated with a specific activation state within the M1-M2 spectrum, identified by specific markers (Murray et al., 2014). Once isolated and seeded on well plates, macrophages can be used as models of the inflammatory response and as phagocytosis assay (Fraser et al., 2009;Westman et al., 2020). As for the in vitro mono-culture of SPCs, well plates are routinely used to culture cells and evaluate properties such as cellular proliferation, differentiation, metabolism and senescence following standard protocols (Groeneveldt et al., 2020). Co-culture systems study the interaction between the immune cells and SPCs to model the inflammatory response in bone healing. Traditional co-culture systems consist of both direct co-culture, where cells are in direct contact with each other on cell culture plastics, and indirect coculture, where transwell inserts are added to culture plates to keep the two cell types separated from each other (Goers et al., 2014). By tuning the pore size of the transwell insert, indirect co-culture models were employed to study the paracrine cell-cell signaling (pore size 0.4 μm, Zhang et al., 2017) or the chemoattractant effect of immune cells on SPCs (pore size 8 μm, Anton et al., 2012). Recent reviews discuss in vitro models of the interaction between SPCs and T cells (Kovach et al., 2015) or macrophages (Maruyama et al., 2020), as well as their implications for bone healing.

Advanced in vitro Models
Besides traditional cell culture plastics, novel in vitro systems enable higher control of the culture environment and cellular interaction. Recent developments of the organ-on-chip technology provide confined engineered microenvironments where biochemical and physical stimuli can be finely tuned over space and time (Zhang et al., 2018). By changing the culture chambers design, organ-onchips can incorporate multiple cell types cultured both in 2D, as the endothelial monolayer (Del Amo et al., 2016), or in 3D, as cells embedded in a hydrogel (Nasello et al., 2020). The optical transparency of organ-on-chip devices facilitates live cell imaging and monitoring. Compared to traditional transwell inserts, organ-on-chips provide a more physiological environment to study the transendothelial migration of inflammatory cells (Han et al., 2012) and the recruitment of SPCs (Eng et al., 2013). In addition, these systems facilitate the application of both chemical (Moreno-Arotzena et al., 2014) and mechanical cues (Middleton et al., 2017) during culture (Occhetta et al., 2019). Therefore, the combination of inflammation-on-chip (Irimia and Wang, 2018) and bone-on-chip (Nasello et al., 2021) would offer a unique alternative to validate in silico models of bone healing by replicating key cellular and environmental interactions of the inflammatory phase.

In vivo Models
Despite ethical concerns, validation with animal models is still an essential step for any preclinical study of both the immune system (Wagar et al., 2018) and the bone repair process (Mills and Simpson, 2012;Lammens et al., 2021). Based on the size, in vivo models are generally divided into small and large animals. Their use depends on the biological process under investigation and the translational stage of the study ( Figure 3A). Here, we discuss the most common small and large animal models used when focusing on the role of the immune system in bone fracture healing. The main results are commented from the modeler's perspective, in view of creating in silico counterparts of these studies. When modelers retrospectively collect data from animal studies to estimate input parameters, they should be aware of the physiological differences between anatomical regions of the skeleton. Besides differences in developmental origins, structural variations in bone composition and direct changes in the biomechanical environment (Cointry et al., 2016), there are regional specializations in cellular composition and differentiation potential. For example, the differentiation potential of skeletal progenitor cells varies between different anatomical sites both in small and large animals (Groeneveldt et al., 2020;Sivaraj et al., 2021).

Small Animal Models
Murine models are widely used to study human diseases and physiology. Despite the differences in the immune system (Mestas and Hughes, 2004) and in the fracture healing process (Haffner-Luntzer et al., 2016) of rodents and humans, murine models can provide clinically relevant results. For instance, murine models were used to validate the clinical observation that fracture healing rate is correlated to higher levels of CD8 + T cells in the peripheral blood (Reinke et al., 2013). By depleting or introducing CD8 + T cells in a mouse model, the authors observed that fracture regeneration was enhanced or impaired, respectively. Therefore, when modeling the immune effect on bone healing, the levels of CD8 + T cells in peripheral blood might be used as a marker of the patient-specific immune reactivity (Reinke et al., 2013).
As for the murine model choice, conventional mouse inbred strains are commonly used since animals share an almost identical genotype, thus leading to higher consistency in the experimental results (Wagar et al., 2018). The lack of genetic variability between individuals is the reason why researchers prefer inbred strains to investigate the fundamental effects of the inflammatory response during fracture healing. The depletion of specific immune cell types, such as macrophages (Schlundt et al., 2018) and T cells (Reinke et al., 2013), was assessed in the mouse inbred strain named C57BL/6N. Moreover, the same inbred strain was used to demonstrate that T and B cells invade the fracture site during the inflammatory phase and the callus mineralization (Könnecke et al., 2014).

Large Animal Models
Large animal models are the most realistic experimental models of human biology and therefore an essential pre-clinical step in translational research (Ribitsch et al., 2020). While nonhuman primates are the most representative model of the human immune system (Wagar et al., 2018), pigs and sheep are normally used to model bone repair since their bone anatomy, mineral composition, regeneration capacity and biomechanical properties are relatively similar to human's (Sparks et al., 2020). Moreover, compared to mice, their immune system is closer to the human one (Lüthje et al., 2018). As a consequence, pigs and sheep are the first choices as large animal models of the inflammatory response in bone fracture healing.
Compared to small animal models, the biological responses of large animal models are more heterogeneous. While small animal models mostly provide mechanistic insights, such as the effect of depleting a specific cell type, research using large animal models tends to explore the complete biological response and the effects on the entire organism, namely the systemic effects. An in vivo study on pigs showed temporal differences in the upregulation of pro-inflammatory cytokines at the fracture site and in the peripheral blood (Horst et al., 2015). Therefore, the validation of in silico models using the cytokine levels in blood as input is intrinsically related to a large animal study.
Another key advantage of using large animal models is the possibility to apply clinically relevant mechanical loads to the fracture site. Schmidt-Bleek et al. (2012) showed that mechanical loads delaying bone healing corresponded to a higher presence of T cells in the fracture site, a prolonged inflammatory signaling in the periosteum and reduced angiogenesis. Hence, large animal models should be chosen to assess the interplay between the immune system, bone repair and mechanical loads.

Laboratory Techniques for Experimental Evaluation
Experimental cell-scale techniques can validate in silico models describing cellular functions and their response to external stimuli. Therefore, this subsection discusses the quantification of molecular mechanisms behind cell processes which could be applied to both in vitro and in vivo experiments. Additionally, live imaging techniques are discussed to calibrate cell invasion parameters with in vitro migration assays. As for extracting tissue-level information, standard imaging methods consist of micro-computer tomography (micro-CT), histology and immunohistochemistry. While their use in preclinical models of bone defects has been recently described elsewhere (Sparks et al., 2020), the present subsection shows examples of tissue-level data extracted from images that could validate in silico models.

Cell-Scale Techniques
The structural and functional characterization of biological molecules belongs to the scientific fields named omics. To explore cellular processes in bone biology, omics technologies characterize, among others, DNA modifications (epigenomics), RNA transcriptions (transcriptomics), protein synthesis (proteomics) and metabolic activity (metabolomics) (Reppe et al., 2017). Transcriptomics, proteomics and metabolomics provide a direct measure of cell survival, proliferation, differentiation and phenotype (Calciolari and Donos, 2020). Therefore, omics technologies can validate in silico models of bone healing by coupling cellular function to tissue adaptation ( Figure 3B). Regarding transcriptomics, RNA sequencing (RNA-seq) technologies measure whole transcriptomes, thus they simultaneously analyze the gene expression profile of thousands of genes (Stark et al., 2019;Calciolari and Donos, 2020). When applied to fracture healing, RNA-seq revealed differences in gene expression associated to skeletal and vascular formation between two mice strains, which was correlated to differences in endochondral bone formation (Grimes et al., 2011). Full and stress fractures revealed different transcriptional profiles during repair, with higher expression of inflammatory and immunerelated genes in full fractures (Coates et al., 2019). In addition, RNA-seq can measure the dynamic changes in the transcriptome. For example, RNA-seq showed that bone marrow stromal cells upregulate pro-inflammatory gene expression during aging, supporting the hypothesis of a regulatory effect on hematopoietic stem cells (Helbling et al., 2019).
As for proteomics, multiplex immunoassays use specific detection antibodies to measure the level of target proteins. Therefore, immunoassays can quantify a large set of inflammatory cytokines from serum or hematoma samples (Horst et al., 2015), as well as cytokines and osteogenic factors synthesized by macrophages and osteoprogenitor cells in vitro (Zhang et al., 2017). Another analytical tool in proteomics is mass spectrometry, which can be used to evaluate the protein composition of the SPC secretome and its variation in a proinflammatory environment (Maffioli et al., 2017).
Concerning metabolomics, cell metabolism is continuously altered in bone healing (Loeffler et al., 2018) and experimental techniques can measure both the metabolites produced and the metabolic pathway activities. Glucose and lactate levels are routinely measured in cell culture media, while glucose consumption and lactate secretion can be calculated by comparison with the unspent medium. An increase in glucose uptake and lactate secretion is associated with the M1 macrophages (Galván-Peña and O'Neill, 2014), but their direct calculation from the cell culture media cannot be related to the specific metabolic pathway producing lactate from glucose. Variations in metabolite levels in the media might be related to faster uptake or lower secretion, rather than to a switch in the metabolic routes. To quantify the activity of each metabolic pathway, glucose is labeled with isotope tracers and incorporated radioactivity is measured. Therefore, isotope tracing reveals differences in glucose uptake for M1 and M2 macrophages (Vats et al., 2006).
It is important to mention that traditional omics technologies consist of analyzing the bulk sample, meaning that the quantified data refers to the whole cell population or tissue. Therefore, bulk omics technologies lose the information regarding RNA transcription, protein synthesis and metabolic activity of individual cells. To maintain the biological information of individual cells, novel advances in the omics field separate and analyze single cells from the population (Barh and Azevedo, 2019). For example, high-throughput techniques for RNA-seq allow to measure the whole transcriptome of single cells (scRNAseq) (Goodwin et al., 2016). By using multiplexed and parallel detection systems, scRNA-seq is generating data that can be used to construct cell atlases from animal and human tissues, both from pathological and physiological conditions (Camp et al., 2018). Although the majority of scRNA-seq methods do not preserve the spatial information of transcriptomic data, novel methods are arising to first localize cells in tissue sections and then sequence RNA (Camp et al., 2018). It is clear that the calibration and validation of in silico models, especially discrete models, would benefit from cell atlases reporting the variations of the transcriptome in bone fractures over space and time.
Besides omics technologies to quantify the molecular mechanisms, migration assays are relevant experimental techniques to measure cell-scale parameters in bone healing. For instance, cell culture inserts can assess SPC migration in vitro and already quantified higher migration capacity under inflammatory conditions (Anton et al., 2012). However, the in vitro environment of a cell culture insert does not represent the 3D extracellular matrix in which cells are embedded in vivo. In search of a more representative structural and biological environment, organ-on-chip systems offer confined 3D culture chambers to monitor the migration of inflammatory and osteoprogenitor cells throughout the experiment (Del Amo et al., 2018;Irimia and Wang, 2018). By tuning the microstructural properties in the culture chamber, such as using fibrin or collagen hydrogels, the organ-on-chip can mimic the different extracellular matrices during bone repair. By changing the cell types and the mechano-chemical stimuli in the device, the organon-chip can selectively identify the role of different factors in cell migration. Therefore, in silico modelers can use cell migration assays in organ-on-chips both to calibrate and validate their predictions.

Tissue-Scale Techniques
Experimental methods to evaluate bone repair at tissue level are widely applied to animal studies and they mostly rely on imaging techniques. Imaging techniques provide qualitative and quantitative information about the analyzed tissue, thus they can validate bone healing in silico models ( Figure 3B). For example, micro-CT provides metric and non-metric parameters of the bone tissue, such as the mineral density of the bone matrix and trabecular morphology (Müller, 2009). Micro-CT imaging has been applied recently for the noninvasive monitoring of fracture healing in mice. By registering time-lapsed scans of the fracture, micro-CT facilitates the assessment of bone parameters throughout the healing process, without altering the callus properties (Wehrle et al., 2019). This imaging technique has already been coupled to FE models of the mechanical in vivo environment, as it substantially contributed to the creation of a personalized bone regeneration model (Tourolle né Betts et al., 2020). Moreover, the combination of micro-CT images and FE models reveals the influence of mechanics on the processes of bone formation and resorption (Birkhold et al., 2014), which can be used to validate in silico models of bone adaptation (Schulte et al., 2013).
While micro-CT imaging allows to quantify the newly formed mineralized tissue, histological sections provide histomorphological parameters of the regenerated bone. The histomorphometrical analysis of the Movat Pentachrome staining quantifies the relative area of bone marrow and connective, cartilaginous and osseous tissue (Schlundt et al., 2018). Additionally, immunohistochemical analyses can be used to stain specific cells in a tissue and quantify their density, thus they can identify the different cell types within the fracture site. The output of immunohistochemical analyses is the fraction of the target cell type, such as CD8 + T cells, M1 macrophages or osteoblasts, in the stained section (Wendler et al., 2019).

TOWARDS THE NEXT GENERATION OF BONE HEALING IN SILICO MODELS
Considering all the computational models of bone healing that have been developed in the last years, it is surprising that almost none of them describe the early inflammatory phase. As the initiator of the bone healing process, inflammation has a considerable impact in the later stages of bone repair: if the inflammatory response is too strongly down-or upregulated, the fracture can result in non-union. To our knowledge, only one in silico model of the inflammatory response in bone healing was developed, which captured the effect of pro-and antiinflammatory pathways on the healing outcome (Kojouharov et al., 2017;Trejo et al., 2019). Although this model laid a strong foundation within the field of computational bone fracture healing, it still has limitations, among them the lack of spatial distribution of cells and cytokines within the healing region. Therefore, this review aims to guide the design and the validation of the next generation of bone fracture healing in silico models, which will include the inflammatory phase.
The biological problem was initially defined by exploring the process of bone healing, with particular attention to the inflammatory phase and its cellular and subcellular components. The inflammatory reaction after bone fracture is a highly complex process, as there is an interplay between different levels (tissue, cellular and subcellular) and systems (musculoskeletal and immune). However, while the mechanobiological activities at cellular and subcellular levels are usually challenging to investigate experimentally, the in silico approach can be employed to unveil the hidden events happening at different time and length scales. Following the current trend of developing hybrid multiscale models (Checa et al., 2011;Carlier et al., 2012;Ceresa et al., 2018;Borgiani et al., 2019) to integrate individual (sub)cellular contributions to tissue dynamics, it seems straightforward that the computational research of the inflammatory phase in bone healing should take advantage from similar methodologies. Hybrid multiscale models benefit from both continuous models, which capture the mechanoregulation of tissue and cellular dynamics, and discrete models, which describe the stochastic interactions at cellular and subcellular level occurring during the immune response.
At the current state-of-the-art, numerous in silico models of the inflammatory response for different organ systems have been developed within the field of computational immunology (see Subsection 3.2). Since the inflammatory response always tends to follow an analogous cascade of events, these in silico models generate the basis to simulate the inflammatory phase in bone healing. Within computational immunology there is a clear preference to use discrete approaches, such as agent-based models or cellular automata, to represent the stochasticity of the immune system. Discrete models can capture better the cascade of cells and subcellular factors that characterize the inflammatory response in bone healing at different time and length scales. However, continuous algorithms capture the dynamics of tissue formation during the subsequent repair and remodeling phases. The development of a comprehensive model that can simulate the mechano-regulation of tissues and the dynamics of large cell populations, while accounting for the probabilistic rules dominating the biological events at subcellular level, requires the combination of continuous and discrete models. For this reason, we believe that the next generation of in silico bone healing models will rely on hybrid approaches to include inflammatory regulation.
In order to guarantee the credibility of in silico models results in a clearly defined context of use, verification, validation and uncertainty quantification analysis (VVUQ) (ASME, 2018; Parvinian et al., 2019) are essential. Verification ensures the accuracy of the model implementation and validation confirms the correspondence between simulation results and experimental reality. The correspondence between computational outputs and physical reality is intrinsically related to in vitro and in vivo experiments; as the in silico modeling of biological processes, like the ones listed in Table 1, requires thorough parameter estimation. Computational modelers should take advantage of different experimental setups able to provide data for the time and length scales simulated. Table 1 presents a proposal to validate certain biological activities happening during the bone healing process using in vitro techniques to replicate specific biological mechanisms in a laboratory, thus providing quantitative data to estimate model parameters. In general, in vitro models and multi-omics approaches can validate in silico models describing signaling pathways involved in cell fate decision or the response of different cell types under external cues. Therefore, in vitro models and omics approaches are highly recommended for the validation of discrete computational models simulating biological events at cellular or subcellular levels. In vivo models and imaging techniques are more suitable to validate continuous or hybrid in silico models describing the biological response at higher scales, such as histomorphometrical parameters or tissue mineralization. As a result, multiscale or hybrid models covering different time and length scales might require both in vitro and in vivo models, as well as both cell and tissue level experimental techniques, for their validation. Furthermore, in vitro experimental studies are performed to calibrate in silico models during their design. The possibility to isolate single biological mechanisms in vitro and introduce them as calibrated parameters within the in silico model allows to simulate behaviors that resemble the ones observed experimentally and investigate their role in the overall outcome of the simulation. The impact of each parameter on the simulation can be quantified with uncertainty quantification methods.
Uncertainty quantification is tested with sensitivity analyses using e.g. Design of Experiments (DOE) or Machine Learning approaches (Mehrian et al., 2018) to assess whether the uncertainty in model assumptions and parameter values does not lead to nonphysiological results. The use of these methodologies to investigate and evaluate the inference of the different parameters can result in valuable information about the most realistic values to describe mechano-biological events. For instance, Isaksson et al. (2008) used DOE to evaluate the significance of multiple factors in bone fracture healing. Parametric uncertainty was addressed by evaluating the outcome of different experiments (simulations, as the study was performed in silico) characterized by organized combinations of parametric values assigned to the factors that describe the bone fracture healing process at cellular level (Isaksson et al., 2008). Another class of optimization techniques is the one used by Steiner et al. (2013), namely evolutionary computation. They calibrated their in silico model by using the Particle Swarm Optimization (PSO) method to achieve the optimal characterization of the mechanical properties of the tissues in a bone fracture healing scenario. The PSO algorithm evaluated combinations of parameters equally distributed in a stochastic way to find the best combination to describe the tissue mechanics (Steiner et al., 2013). Machine Learning techniques can also be used to evaluate the best value fitting of specific parameters or to categorize certain outputs. An Artificial Neural Network was used by Cilla et al. (2017) to evaluate the geometrical features for the design of a patient-specific short-stem hip implant to contrast the mechanical side effect of prosthetic stress-shielding (Cilla et al., 2017).
The inflammatory response in bone fracture healing has a noteworthy complexity, but in silico models help us to understand the principles regulating the diverse events occurring at tissue, cellular and subcellular level. Certainly, the experimental validation of such in silico models is mandatory if we aim to go from bench to bedside. With this review, we aimed to highlight the potential of using multiscale in silico approaches to tackle bone healing intricacy. Based on the current state-of-the-art, we conclude that hybrid models are particularly suited to simulate adequately the multiscale course of events of the inflammatory phase and its overall role in the healing outcome. We furthermore described possible in vitro and in vivo methodologies that can be employed to experimentally calibrate the parametric description of the in silico model during its development and, afterwards, to validate the computational results and support their bench to bed transition. We believe that the next generation of in silico models of bone regeneration should account for inflammatory events to guarantee a more realistic investigation of the process, favoring its employment within a clinical context.

AUTHOR CONTRIBUTIONS
LL-G, EB, and GN drafted the manuscript. All authors read and revised the manuscript and approved its final content.