In this study, a coupled compu tational modelling framework for bone fracture repair is presented that enables predictions of both healing and remod- elling phases of the fracture region and is used to investigate the role of an inter- nal fixation plate on the long -term healing performance of a fracture tibia under a range of different conditions. It was found that introduction of a titanium plate allowed the tibia to undergo successful healing at higher loading conditions and fracture gaps, compared to the non -plated versions. While these plated cases showed faster rates of repair in the healing phase, their performance was substan- tially different once they entered the remodelling phase, with substantial regions of stress shielding predicted. This framework is one of the few implementations of both fracture healing and remodelling phases of bone repair and includes sev- eral innovative approaches to smoothing, time -averaging and time incrementa- tion in its implementation, thereby avoiding any unwanted abrupt changes be- tween tissue phenotypes. T his provides a better representation of tissue develop- ment in the fracture site when compared to fracture healing models alone and provides a suitable platform to investigate the long -term performance of ortho- paedic fixation devices. This would enable the more effective design of perma- nent fixation devices and optimisation of the spatial and temporal performance of bioabsorbable implants. Keywords: Bone Fracture Healing; Finite Element Modelling; Orthopaedic Im- plants 1. Introduction Secondary bone fracture repair is a four -stage process [1] that begins with an inflammatory response and haematoma formation, resulting in granulation tissue formation 3 -7 days after the initial fracture. After 2 -4 weeks, a fibro -cartilaginous soft callus is produced from the g ranulation tissue, which subsequently leads to the development of an ossified hard callus of woven bone surrounding the fracture region after 2 -4 months. Over the next several months to years, remodelling of the woven bone takes place to achieve the origin al pre -fracture morphology [1], [2] . Assuming the fracture site receives an adequate supply of blood, hormones and growth factors [3], fracture healing is largely dictated by Copyright: © 2021 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/license s/by/4.0/). 2 the biomechanical conditions in the fracture callus [4], [5] . To enhance th e repair process, various types of internal fixation devices are available, whose role is to stabilise the fracture and minimize inter - fragmentary movement in the early phases of healing. Many complications in bone repair, such as non - union, implant failur e, and delayed healing, are derived from instability in the fracture region [6]. However, fixation devices themselves can cause late -stage problems in the repair process, as bone resorption can take place around overly stiff fixators in the remodelling p hase, which can increase the likelihood of implant failure [7]. The ideal fixation device should reduce inter -fragmentary movement in the early phases of healing to minimise damage to the callus tissue but should enable load -transfer to the bone itself o nce the remodelling phase has initiated to avoid unwanted effects of stress shielding [8]. Understanding the biomechanical performance of these devices represents a challenging mechano - biological problem that needs to consider in detail the device -tissue interaction throughout all stages of the bone repair process. Computational models have been used extensively to understand bone fracture healing, whereby the iterative finite element method has been used to determine the mechanical environment at the si te of fracture and local tissue formation proceeds according to appropriate mechano -regulation theories [1], [2], [9], [10] . Typically, these models account for a range of biophysical stimuli at the fracture site, including combinations of strain, stress , fluid flow, streaming potentials and acceleration [11], that drive biological processes and cellular activities, including migration and proliferation of mesenchymal stem cells (MSCs) from the external environment, their differentiation into fibroblast s, chondrocytes or osteoblasts and tissue formation. While these models require careful calibration of numerical parameters, they have provided a detailed understanding of how local factors, such as mechanical stimuli and cellular activities, affect tissue production following a bone fracture without fixation [1], [2], [9], [10] . Several studies have also used similar frameworks to explore the role of implant design and fixation conditions on fracture healing response [16], with studies that identify re ductions in healing time are possible with fixators of increasing stiffness [16]. However, a limitation of such modelling approaches is that the final remodelling phase of bone repair is rarely fully considered. Instead, these models tend to examine the role of fixation on inter -fragmentary movement in the early phases of healing, but do not consider wider remodelling processes that take place in both the callus and the adjacent cortex that take place following a successful union over longer time periods [1], [16]– [18] Healthy bone tissue actively remodels its apparent bone density in response to the local mechanical environment. The constant remodelling of bone from the cellular to organ scales means that bone density is optimised as a function of the loading environment according to Wolff’s law. Bone remodelling algorithms are widely used to provide numerical estimates of bone’s apparent density due to daily loading according to the theory of mechanostat [19]. Numerous bone remodelling theories have been proposed, with theories using strain -derived [20]–[22] and/or damage -based [23]–[29] stimuli to describe cellular activity [30], [31] and drive the remodelling process. While these theories support the concept that bone remodelling is a targete d mechanism to simultaneously preserve bone mass and prevent fractures, bone remodelling also plays a vital role in the final stage of the fracture repair process, as woven bone is gradually replaced with orientated lamellar bone to return to original pre -fracture morphology. However, these remodelling processes have rarely been considered in existing healing 3 algorithms, which means that the longer -term tissue responses have not been fully explored. The ability to predict the response of bone tissue to adap tations in mechanical loading [32], [33] is of vital significance when considering how the introduction of load -bearing implant will affect the structural architecture. Identifying a fixation device that has optimal stability in the early phases fracture repair and an appropriate design to minimise stress shielding requires computational models of fracture repair that consider both fracture healing and subsequent remodelling responses. These approaches could be particularly relevant given the ongoing deve lopment of bioabsorbable fixation devices, which provide the potential to control both spatial and temporal load -bearing performance during healing. The objective of this study is to develop a coupled modelling framework for bone fracture repair that pred icts both fracture healing and remodelling phases of the tissue over a period of 12 months. During the fracture healing phase, the model predicts tissue production based on a combination of the biophysical stimuli and cellular concentration in the fracture callus. The remodelling phase is only initiated once a suitable fracture union is achieved and considers both damage and strain -energy density (SED) based tissue remodelling. The coupled modelling framework for bone fracture repair is used to investigate the role of an internal fixation plate on the long -term healing performance under a range of loading conditions. 2. Materials and Methods 2.1 Model Geometry Three -dimensional finite element models of non -plated and plated long -bone tibia fractures were developed, according to dimensions provided in Figure 1. The tibia was simplified as a smooth co -axial cylinder composed of cortical bone and bone marrow at the centre (see Figure 1a), with the shape and dimensions based on previous studies [9], [34], [35 ] . The fracture callus was divided into three zones, (1) outer callus, (2) callus focus and (3) inner callus . Several fracture gaps ( x=0.75mm, 2.5mm, 5mm, 8mm and 12mm) were considered. For the plated models, six self -locking screws fastened a titanium bone plate to the tibia [36] (see Figure 1c and d). The cortical bone, callus, bone marrow and screw -plate system were modelled with 32,076, 17,750 and 15,710 8 -noded linear brick, trilinear displacement, pore pressure and temperature, reduced integration elements (C3D8RPT). Figure 1: Geometric information on human tibia implanted and external fixator; (A) bone plate; (B) callus re- gions; (C) top -down view of plate without pins; (D) side view of plate with pins. All measurements were rec- orded in millimetr es. 4 2.2 Bone Repair Model 2.2.1 General Overview The bone repair model is represented schematically in Figure 2, whereby an iterative approach that considers both (a) fracture healing and (b) remodelling phases of bone repair was developed. For the fractu re healing component (see Figure 2a), a fracture callus comprised of granulation tissue is infiltrated by considering angiogenesis and migration of mesenchymal stem cells (MSCs) from the adjacent undamaged cortical ends (see Section 2.2.2). As the patient gradually loads the fractured bone, biophysical stimuli acting on the cells drive differentiation and tissue production process, with gradual changes from one phenotype to another enabled by a smoothing process within the fracture healing algorithm (See Se ction 2.2.3). Once adequate woven bone has formed in the “callus focus” region (see Figure 2b), the entire fracture callus automatically enters the bone remodelling algorithm, whereby both damage and strain -energy density (SED) based remodelling is predict ed. The remodelling algorithm updates the bone material properties according to one of four processes; (1) damage -based remodelling; (2) positive SED remodelling; (3) negative SED remodelling; (4) no remodelling (See Section 2.2.4). The apparent bone densi ty continuously remodels to achieve the homeostatic bone energy density. The model is implemented in the Abaqus/Standard finite element framework through a series of user -defined field Subroutines (USDFLD). Figure 2: Iterative calculations for tissue development for bone fracture repair and bone remodelling; (a) Bone fracture repair algorithm and (b) SED/damage -based bone remodelling. 2.2.2 Blood Vessel Growth/Cell Migration and Proliferation A diffusion model was implemented in the fracture callus t o model the infiltration of MSCs and the development of blood vessels according to, 5 𝑑𝑐𝑚 𝑑𝑡=𝐷𝛥2𝑐𝑚 [1] where C is the cell concentration of MSCs, D is the mass diffusion coefficient. This diffusion model was implemented in Abaqus using a thermal analogy, given the similarity between the therma l diffusion and mass diffusion (Fick's 2nd law) equation. This approach has been applied by other studies, with further details of the implementation process available in [2]. Briefly, the diffusion model assumes that the movement of the cells is random and non -directional, with the local relative cellular concentration ranging from 0≤ 𝑐𝑐𝑒𝑙𝑙≤1. The source of the progenitor cells was the surrounding soft tissue, bone marrow and inner cambial of the periosteum (see Figure 3A) [37]. At the beginning of simulations ( 𝑡= 0), the callus does not contain any MSCs (𝑐𝑚=0) and cons ists of only granulation tissue. The cell concentration 𝑐𝑚 changes with respect to time according to Equation 1. Simulations were initially ran with a range of diffusion coefficients ( D), with a diffusion coefficient of D=0.3mm3/day yielding the best results for cellular confluence following 112 days [35]. Cell concentration at the integration points are calculated through interpolating nodal cell concentrations. Previous studies [2], [38] applied the maximum cellular concentration (𝑐𝑚=1,𝑇𝑛𝑜𝑑𝑎𝑙 =1,@ 𝑡=0) at the boundaries shown in Figure 3A. However, this results in the integration points of these elements having a cellular concentration equal to half of the applied cellular concentration at time zero (𝑐𝑚=0.5,@ 𝑡=0), and causes these elements to stiffen much quicker than the rest of the elements within the fracture callus. Therefore, using Equation 2, a temporal non -linear magnitude was created (Figure 3b) and applied to the fracture callus (Figure 3a). (c m) 𝑐𝑚(𝑖+1)=𝑐𝑏𝑐(𝑖)+𝐷(𝑐𝑚𝑎𝑥 −𝑐𝑚(𝑖))∗(𝑑𝑡) [2] Figure 3: Overview of diffusion model used to describe the infiltration of MSCs/development of blood vessels; (A) MSCs invading the fracture callus from the surrounding soft tissue, bone marrow and inner c ambial of the periosteum; (B) Cellular concentration magnitude applied to surface shown in Fig.3(A) . 2.2.3 Fracture Healing Mechano -Regulation Algorithm The mechano -regulation theory described by Prendergast [15] was used to govern MSC differentiation an d tissue phenotype under mechanical loading within the callus tissue (see Figure 4a). It was assumed that the tissue was biphasic, and completely saturated, such that the sum of the volume fractions of the solid ( 𝑉𝑠) and fluid ( 𝑉𝑓) phase must equal one. 6 𝑉𝑠+𝑉𝑓=1 [3] Mechanical loading produced a deviatoric strain (shear deformation) and a fluid flow throughout the fracture callus. The effective pore fluid velocity ( 𝑤𝑓) was calculated as the velocity of the fluid phase (𝑣𝑓) with respect to the solid phase (𝑣𝑠). 𝑤𝑓= 𝑣𝑓−𝑣𝑠 [4] The deviatoric strain (𝜀𝑑𝑠𝑠) was calculated using Equation 5, which considers this term as a function of the principal strains ( ε1, ε 2, ε 3), while the magnitude of the fluid flow velocity was described in Equation 4. The local elemental stimulus ( S) (Equation 7) was defined as the sum of the equivalent deviatoric strain and the interstitial fluid velocity in which 𝑎=0.0375 and 𝑏=3𝜇𝑚/𝑠 [39]. 𝑤𝑚𝑎𝑔𝑓=√𝑤𝑓.𝑤𝑓 [5] 𝜀𝑑𝑠𝑠= 2 3√(𝜀1−𝜀2)2+(𝜀2−𝜀3)2+(𝜀3−𝜀1)2 [6] 𝑆=(𝜀𝑑𝑠𝑠) 𝑎+𝑤𝑚𝑎𝑔𝑓 𝑏 [7] Figure 4: (A) The mechano -regulation model based on deviatoric strain and fluid velocity. (B) Smoothed mecha- nobiological regulation based on Sapotnick et al. (2015) . The local stimulus within an element controlled whether the MSCs w ithin that element differentiated into fibroblasts, chondrocytes or osteoblasts. As a rule, high levels of mechanical stimuli (S >3) resulted in the formation of fibrous tissue, intermediate levels of stimulation (1< S<3) promoted cartilaginous tissue, lowe r levels of stimulus resulted in formation of immature bone 7 (0.2667 <S<1), while woven bone was formed between when the stimulus was 0.01<S< 0.2667 . For very low stimuli (0.01<S) callus resorption occurs, whereby osteoclasts removed non -loaded bone tissue [2]. These mechano -regulation rules were implemented considering fuzzy zones [40], that allowed for gradual changes in tissue properties as differentiation occurred, whereby two tissue phenotypes could be produced simultaneously, allowing for smoother t ransitions between cell phenotypes. Sapotnick and Nackenhorst altered the range of their fuzzy zones from 50% to 400% with minimal effects on healing outcomes. A relatively conservative approach was taken when determining the width of the fuzzy zones, wher eby the range was determined by taking 10% of the threshold stimulus (S=3,S=1,S=0.2667 ,S=0.01). Tissue membership (𝑇𝑀) was determined by the biophysical stimulus (Figure 4(A)). The tissue membership level -controlled cell differentiation and tissue production. The local TM for fibrous, cartilage, immature bone and mature bone ranged from 0≤ 𝑇𝑀 𝑓,𝑐,𝑖𝑚,𝑚𝑏≤1. The evolution of fibroblasts, chondrocytes, immature bone cells and mature cells and their associated tissues are described by Equations 8 -11 [41]. Here 𝑐𝑚𝑏,𝑐𝑖𝑏,𝑐𝑐,𝑐𝑓 and 𝑐𝑚 are normalized mesenchymal, fibroblast, chondrocyte, immature bone, and mature bone cell densities, respectively, 𝐹𝑚𝑏,𝐹𝑖𝑏,𝐹𝑐 and 𝐹𝑓 are cell differentiation rates. 𝑑𝑐𝑓 𝑑𝑡=(𝐹𝑓(1−𝑐𝑓)𝑐𝑚−𝐹𝑐(1−𝑐𝑐)𝑐𝑚−𝐹𝑖𝑏(1−𝑐𝑖𝑏)𝑐𝑚−𝐹𝑚𝑏(1−𝑐𝑚𝑏)𝑐𝑚)(𝑇𝑀𝑓) [8] 𝑑𝑐𝑐 𝑑𝑡=(𝐹𝑐(1−𝑐𝑐)(𝑐𝑚+𝑐𝑓)−𝐹𝑖𝑏(1−𝑐𝑖𝑏)𝑐𝑐−𝐹𝑚𝑏(1−𝑐𝑚𝑏)𝑐𝑐)(𝑇𝑀𝑐) [9] 𝑑𝑐𝑖𝑚 𝑑𝑡=𝐹𝑖𝑚(1−𝑐𝑖𝑚)(𝑐𝑚+𝑐𝑓 +𝑐𝑐)−𝐹𝑚𝑏(1−𝑐𝑚𝑏)𝑐𝑐(𝑇𝑀𝑖𝑚) [10] 𝑑𝑐𝑚𝑏 𝑑𝑡=𝐹𝑚𝑏(1−𝑐𝑚𝑏)(𝑐𝑚+𝑐𝑓 +𝑐𝑐+𝑐𝑖𝑚) (𝑇𝑀𝑚𝑏) [11] Tissue production and replacement is regulated by the corresponding cells, tissues themselves and TM. The normalized volume fraction of fibrous tissue ( 𝑚𝑓),, cartilage ( 𝑚𝑐),, immature bone ( 𝑚𝑖𝑏), and mature bone ( 𝑚𝑚𝑏) are described by (Equation 12 -15). 𝐷𝑚𝑏,𝐷𝑖𝑏 and 𝐷𝑐 are tissue absorption rates whereby 𝐷𝑚𝑏=𝑄𝑚𝑏,𝐷𝑖𝑚=𝑄𝑖𝑚 and 𝐷𝑐=𝑄𝑐. 𝑚𝑡𝑜𝑡=𝑚𝑓+𝑚𝑐+𝑚𝑖𝑏+𝑚𝑚𝑏 describes the volume fraction of all tissue except granulation tissue, which volume fraction is 𝑚𝑔𝑟𝑎𝑛 =1−𝑚𝑡𝑜𝑡. 𝑑𝑚𝑚𝑏 𝑑𝑡=𝑄𝑚𝑏(1−𝑚𝑚𝑏)(𝑐𝑚𝑏) (𝑃𝑀𝐿𝑚𝑏) [12] 𝑑𝑚𝑖𝑏 𝑑𝑡=𝑄𝑖𝑏(1−𝑚𝑏− 𝑚𝑖𝑏)(𝑐𝑖𝑏)(𝑃𝑀𝐿𝑖𝑏)−𝐷𝑚𝑏𝑐𝑚𝑏𝑚𝑖𝑏𝑚𝑡𝑜𝑡 [13] 𝑑𝑚𝑐 𝑑𝑡=𝑄𝑐(1−𝑚𝑐− 𝑚𝑖𝑏− 𝑚𝑐)(𝑐𝑐)(𝑃𝑀𝐿𝑐)−(𝐷𝑚𝑏𝑐𝑚𝑏−𝐷𝑖𝑏𝑐𝑖𝑏)𝑚𝑐𝑚𝑡𝑜𝑡 [14] 𝑑𝑚𝑓 𝑑𝑡=𝑄𝑓(1−𝑚𝑐− 𝑚𝑖𝑏− 𝑚𝑐− 𝑚𝑓)(𝑐𝑓)(𝑃𝑀𝐿𝑓)−(𝐷𝑚𝑏𝑐𝑚𝑏−𝐷𝑖𝑏𝑐𝑖𝑏−𝐷𝑐𝑐𝑐)𝑚𝑓𝑚𝑡𝑜𝑡 [15] 8 A time control (Δt𝑓𝑟𝑎𝑐𝑡𝑢𝑟𝑒 ) within the fracture healing algori thm was implemented to ensure the change in tissue volume fraction did not exceed the maximum threshold value (Δm𝑚𝑎𝑥 =0.05) for any given time increment . If (𝑑𝑚𝑚𝑏,𝑖𝑏,𝑐,𝑓>Δm𝑚𝑎𝑥) the time increment was updated according to Equation 16. Δt𝑓𝑟𝑎𝑐𝑡𝑢𝑟𝑒 =Δm𝑚𝑎𝑥 (𝑑𝑚𝑚𝑏,𝑖𝑏,𝑐,𝑓), 𝑑𝑡𝑛𝑒𝑤 =Δt𝑓𝑟𝑎𝑐𝑡𝑢𝑟𝑒 ∗𝑑𝑡𝑜𝑙𝑑 [16] The rule of mixtures (RoM) was used to calculate the predicted elastic modulus of the elements within the fracture callus for each incre ment using Equation 17. The RoM calculated the predicted elastic modulus before smoothing occurred, 𝐸𝑓𝑟𝑎𝑐𝑡𝑢𝑟𝑒 =((𝑚𝑔𝑟𝑎𝑛)∗(𝐸𝑔𝑟𝑎𝑛)+(𝑚𝑓)∗(𝐸𝑓)+(𝑚𝑐)∗(𝐸𝑐)+ (𝑚𝑖𝑏)∗(𝐸𝑖𝑏)+(𝑚𝑚𝑏)∗(𝐸𝑚𝑏)) [17] The implementation of this mechano -regulation algorithm also included a further smoothing strategy, whereby the elastic modulus ( 𝐸𝑛+1) was calculated by obtaining a time averaged value of the elastic modulus from the previous 10 days according to, 𝐸𝑡+𝑑𝑡= 1 10∑ (𝐸𝑡)∗(𝑑𝑡𝑡)𝑡=𝑡𝑖 𝑡=𝑡𝑖−10 [18] Many previous studies have implemented time averaging using a fixed increment size [2], [42], [43] . However, using a variable time step with custom time controls allows for greater acc uracy and versatility as the time step (Δt𝑓𝑟𝑎𝑐𝑡𝑢𝑟𝑒 ) is controlled by changes in cellular concentration and prevents abrupt changes in cell/tissue concentrations. Fracture healing continued within the callus until the “callus focus" region was co mpletely connected by differentiated woven bone. Following this, the fracture callus underwent bone remodelling whereby strain energy density (SED)/micro -damage were the controlling stimuli. 2.2.4 Coupled Strain/Damage - based Remodelling The coupled strain - & damage -adaptive algorithm outlined by McNamara and Prendergast [26] was used to simulate bone remodelling. Strain -based remodelling used strain energy density (SED) as the remodelling stimulus ( 𝜓=𝑈/𝜌), where 𝑈 is the strain energy, 𝜌 is the apparent bone density. A site - specific approach was implemented whereby each bone element tried to reach a homeostatic state of SED. The homeostatic reference point for SED was calculated according to, 𝜓𝑈=𝜎𝑈2 (2)(𝐸𝑈)(𝜌𝑈) [19] Where 𝜎𝑈,𝐸𝑈, 𝜌𝑈 are the homeostatic values for stress, elastic modulus, and density, which have values of 6.6MPa, 17000MPa and 1.649g/cm3 respectively. The homeostatic SED reference point Eq.19 is typically calculated at 𝜓𝑅𝑒𝑓 (Figure 5(A)). This however creates can cause difficulties as it creates a dependency between the final converged solution and the initial state and width of lazy zone (𝑤) [44]. Bone elements remodel towards the u pper (𝜓𝑈) and lower (𝜓𝐿) SED values, therefore by calculating the homeostatic SED at 𝜓𝑈, the model will remodel towards homeostatic elastic modulus under homeostatic loading conditions. The lower SED threshold was calculated using (Equation (2 0)). 𝜓𝐿=(𝜓𝑈)∗(1−𝑤) [20] 9 The width (𝑤) value controlled the size of the lazy zone, and the degree of reabsorption. The width (𝑤) had no effect on the value of 𝜓𝑈 as it has no dependency. In line with previous studies, a width of 0.75 was chosen [45]. Bone remodelling was simulated as a change in apparent bone density. Carter and Hayes found that the elastic modulus of cortical and trabecular bone is approximately proportional to the cube and square respectively of the apparent bone den sity [46], this relationship between density (g/cm3) and elastic modulus (MPa) is described by, 𝐸=3790 𝜌3 [21] The rate of remodelling is affected by the amount of free surface area, as remodelling occurs at the pore surfaces. Martin [47] determined the internal free surface area as a function of the apparent density 𝐴(𝜌). The internal free surface area per unit volume 𝑎(𝜌)=𝐴(𝜌)/𝑉 is shown in Figure 5B. At maximum bone density 𝜌𝑚𝑎𝑥=1.73𝑔/𝑐𝑚3 and the minimum bone density 𝜌𝑚𝑖𝑛=0.01𝑔/𝑐𝑚3, no bone remodelling occurs, as the free surface area equals zero. The remodel ling rules and formulae that determine the rate of remodelling are shown in Figure 5C. Figure 5: (A) Bone remodelling algorithm using SED as the mechanical stimulus for the implant;(B) Surface area density of bone (𝑎(𝜌)) as a function of bone density [48]. A strain remodelling rate was determined from 𝐶𝑠=𝑎(𝜌)𝜏, where 𝜏 is a time constant with a value of 𝜏=130 𝑔2𝑚𝑚−2𝐽−1 per month, and 𝑎(𝜌) is the specific surface area [49]. Time controlled incrementation (𝛥𝑡𝑆𝐸𝐷) was implemented such that for any given iteration the maximum allowable density change (Δρ𝑚𝑎𝑥) was 0.175 𝑔/𝑐𝑚3 per month . If the density change for any given iteration exceeded the allowable maximum, then time increment (𝑑𝑡𝑛𝑒𝑤) was reduced according to, 𝛥𝑡𝑆𝐸𝐷 = Δρ𝑚𝑎𝑥 (𝑎(ρ)∗τ∗(S−𝑆𝑈,𝐿)∗𝑑𝑡), 𝑑𝑡𝑛𝑒𝑤 =𝛥𝑡𝑆𝐸𝐷∗𝑑𝑡𝑜𝑙𝑑 [22] Damaged bone tissue within bone is constantly removed and replaced with healthy bone tissue. At equilibrium, damage (𝜔𝑅𝐸) exists within bone due to the cyclic loading of bone tissue. Damage is calculated within each element using the remaining life approach, whereby undamaged bone is denoted by 𝜔=0, while fracture bone is denoted by 𝜔=1 [50]. It is assumed within healthy bone tissue that damaged is removed at a constant rate which is denoted by 𝜔̇𝑅𝐸. If the damage formation 10 rate 𝜔̇ exceeds the repair rate 𝜔̇𝑅𝐸 damage is accumulated within the bone tissue it will create damage remodelling stimulus 𝛥𝜔 as shown in Eq. 18. 𝛥𝜔= ∫(𝑡 0𝜔̇− 𝜔̇𝑅𝐸)𝑑𝑡 [23] The damage formation rate was assumed to be linear and was implemented using Miner’s rule [51], which determined the damage formation rate 𝜔̇ and repair rate 𝜔̇𝑅𝐸. Previous experimental dat a by Carter et al. [52] provides the number of cycles to failure (𝑁𝑓) for several different stress levels. Using Equations 24 & 25, the damage and repair rates were calculated, where 𝜎 is the stress in MPa, T the temperature (37 °c) and H, J, T and M are empirical constants. 𝜔̇= 1 𝑁𝑓 [24] log(𝑁𝑓)=𝐻𝑙𝑜𝑔 (𝜎)+𝐽𝑇+𝐾𝜌 +𝑀 [25] Once the accumulated damage exceeds a critical damage threshold ( 𝜔𝑐𝑟𝑖𝑡), damage -based resorption occurred. The critical damage level was determine d at 3500 μɛ, with 𝜓𝑈 located at 1500 μɛ. Therefore, the critical damage threshold is situated 2000 μɛ above 𝜓𝑢 or approximately 1.3 times the strain of 𝜓𝑈. Damage can exist within the bone without activating damage -based reabsorption, until the damage critical 𝜔𝑐𝑟𝑖𝑡 damage threshold is exceeded. Similarly, to the SED -based remodelling a damage -based time control ( 𝛥𝑡𝑑𝑎𝑚𝑎𝑔𝑒 ) was implemented (Equation 26). A damage rate constant of 𝐶𝑑=0.1𝑒8 was used [49]. Simulatio ns considered cases with damage ( D) and without damage (ND) based remodelling for the plated models 𝛥𝑡𝑑𝑎𝑚𝑎𝑔𝑒 = 0.00625 (𝑎(ρ)∗𝐶𝑑∗N∗𝜔), 𝑑𝑡𝑛𝑒𝑤 =𝛥𝑡𝑑𝑎𝑚𝑎𝑔𝑒 ∗𝑑𝑡𝑜𝑙𝑑 [26] Within combined fracture healing/remodelling al gorithm the smallest value of 𝛥𝑡𝑓𝑟𝑎𝑐𝑡𝑢𝑟𝑒 ,𝛥𝑡𝑆𝐸𝐷 𝑎𝑛𝑑 𝛥𝑡𝑑𝑎𝑚𝑎𝑔𝑒 was always implemented in order to determine the new time increment 𝑑𝑡𝑛𝑒𝑤. 2.5 Boundary Conditions Post -operative loading of the tibia contained four di fferent loading phases (1) indirect loading, (2) walking with crutches, (3) imperfect gait and (4) perfect gait (See Figure 6) [43]. The body weight of a healthy patient was assumed to be 70kg. An initial load of 10% of body weight (BW) was applied for the first 28 days in order to simulate indirect loading. From days 28 to 56 this the tibia was loaded to simulate the patient walking with crutches, gradually transferring more load through the injured tibia. The load transferred through the tibia continual ly increased between days 56 -112 as the patient walked with an imperfect gait. Following 112 days, a normal gait of 300% BW was implemented. The initial loading conditions (LC) were altered between 0 – 112 days, according to several different conditions described by LC(A) – LC(H), as shown in Figure 6. The distal end of the tibia was fixed in all directions [43]. 11 Figure 6: Post -operative loading conditions. 2.6 Material Properties During the healing process, the MSCs within the fracture callus differen tiate into fibroblasts, chondrocytes and osteoblasts [53]. All materials are assumed to be isotropic. The material properties of the different tissues used in the bone fracture repair and bone remodelling algorithm are listed in Table 1. Table 1. – Mate rial properties for bone plate and bone. Elastic Modulus [MPa] Permeability [m4/Ns] Poisson’s Ratio Initial Porosity Cortical bone 17000 [10] 1E-17 0.37 0.04 Marrow 2[38] 1E-14 0.167 0.8 Granulation tissue 0.2 [10] 1E-14 0.167 0.8 Fibrous tis sue 0.2 -5 [10] 1E-14 0.167 0.8 Cartilage 5 – 500 [10] 5E-15 0.167 0.8 Immature bone 500 – 1,000 [10] 1E-13 0.3 0.8 Intermediate bone 1,000 – 2,000 [10] 3.7E -13 0.3 0.8 Mature bone 2,000 -6,000 [10] 3.7E -13 0.3 0.8 Titanium (Grade 2) 110,000 - - - 2.7 Parameter Study A 3mm mid -diaphyseal osteotomy was induced in the tibia of a sheep and stabilised with a unilateral external fixator. An initial displacement of 1.13mm was induced through the application of a 500N force. Calibration of the model produced the following model parameters: cell differentiation rates 𝐹𝑚𝑏 = 0.15day-1, 𝐹𝑖𝑏 = 0.15day-1, 𝐹𝑐 = 0.3day-1, 𝐹𝑓 = 0.02day-1, tissue production rates 𝑄𝑚𝑏 = 0.1day-1, 𝑄𝑖𝑏 = 0.1day-1, 𝑄𝑐 = 0.2day-1 and 𝑄𝑓 = 0.12day-1. Using these parameters, the model was accurately able to reproduce the interfragmentary strain obtained by Claes et al. [3]. Here, the model does not predict the formation of a bony bridge in the osteotomy gap, leaving a portion of the callus filled with fibrous and 12 cartilage tissue. Mesh sensitivity for fracture healing was carried out using a global seed ( GS) of 1.5, 1 and 0.5, which resulted in 4256, 12547 and 91256 elements. Element size had minimal on the healing outcome as shown i n Figure 7.A. For all simulations, a global seed of 1.5 was implemented. Figure 7: (A) Predicted IFM by FE -model versus healing time of an osteotomy of a sheep. (B) Histology sec- tion of a well -healed mouse at day 26. (C) Histology section of a mouse w ith mixed healing at day 26. Stained with Hall’s and Brunt’s Quadruple (HBQ) stain and false -coloured. Blue = cartilage, yellow = trabecular bone, purple = fibrous/amorphous tissue, red = cortical bone, black/white = bone marrow [54]. 3. Results 3.1 Fracture Repair - Effect of Loading Conditions (Non -Plated) Figure 8A and B shows the predicted mean elastic modulus in the callus focus region (e.g. between both cortex ends) over the course of 365 days, based on the coupled fracture healing and bone remodelling algorithms, under a range of different loading conditions. For the non -plated 3mm fracture gap, the rate of healing was reduced when higher loading regimes were applied. Here, i t is shown that the formation of a bony bridge occurred at 16.6 days for LC(A), which was the smallest loading regime, resulting in the fastest rate of healing with fracture union achieved at 120 days. For increasing loading magnitudes, the time to form a bony bridge and fracture union was progressively longer. For excessive loading conditions in LC(D) and LC(E), there is little or no increase in elastic modulus over the first 120 days. The higher loads inhibited the formation of a bony bridge, preventing t he stabilization of fracture callus region and resulted in the formation of a non -union. The coupled fracture repair framework only allows the callus to enter SED based remodelling once a fracture union had formed between the cortical ends. During remodell ing, the SED/microdamage based framework caused the outer and inner callus to be predominately negatively remodelled and reabsorbed. Meanwhile, the callus focus region underwent positive remodelling as shown in Figure 8(E). The undamaged cortical bone init ially negatively remodelled under the reduced physiological load, however as loading returns to normal levels, the undamaged cortical bone was positively remodelled and returned to homeostatic levels (see Figure 8(E)). 13 Figure 8 – Elastic modulus of fractured tibia over a 1 -year time period with a constant fracture gap and different loading conditions; (A) fracture callus healing performance over 1 year; (B) fracture callus healing pe rformance over 112 days; (C – E) contour plots of fracture callus; (F) SED/micro -damage based remodelling. 3.2 Fracture Repair & Remodelling - Effect of Fracture Gap Figure 9 shows the predicted elastic modulus in the callus over the 365 days, based on the coupled fracture repair framework for loading condition LC(A) and fracture gaps of 3mm, 4mm, 6mm and 8mm. Here, the fracture healing performance was directly related to the fracture gap size, with the 3mm showing the fastest rate of healing. Increasing th e fracture gap size slowed the rate of fracture healing as there was a slower infiltration of cells to fracture site, in addition to an increased biophysical stimulus due to geometry. Remodelling proceeded like before, with the inner and outer callus under going negative remodelling, while the callus focus region underwent positive remodelling. Where a union is successfully formed, the elastic modulus of the remodelled tissue is very similar across all cases after 365 days. 14 Figure 9: Elastic modulus of fractured tibia over a 1 -year time period with a constant loading condition and dif- ferent fracture gaps; (A) fracture callus healing performance over 1 year; (B) fracture callus healing performance over 112 days; (C – E) contour plots of fracture callus; (F) SED/micro -damage based remodelling. 3.3 Fracture Repair – Effect of Plating The effects of the implantation of a titanium plate on fracture healing and remodelling was investigated for LC(A), LC(F) –LC(H). Figure 8B shows that the predicted elastic modulus in the callus during the 365 days of repair for these loading conditions, while Figure 9B shows the predicted elastic modulus during repair for differe nt fracture gap sizes for LC(A). Figure 8B shows that the introduction of a titanium plate allowed the tibia to undergo successful healing at much higher loading conditions, compared to the un -plated versions. Non -union occurred for the un -plated case unde r LC(C) (0.33MPa), while the plated version successfully healed up to LC(H)(6.6MPa), which had substantially higher magnitudes applied (see Figure 6). The introduction of a stiff plated increased the overall rate of fracture repair, which is clear when com paring the plated/un -plated models utilizing LC(A). Figure 10 shows the healing profile of the callus itself, where the presence of the plate affects the spatial 15 progression of healing. Tissues positioned closer to the plate stiffened and was reabsorbed qu icker than tissue positioned further away for the plate. While the plated cases showed faster healing rates, their performance was substantially different once they entered the remodelling phase, with all plated cases remodelling to a much lower elastic mo dulus than any of the un -plated loading cases (see Figure 8). The introduction of a titanium plate allowed for successful fracture healing to occur at much larger fracture gaps, although with a delayed healing performance with increasing fracture gaps (Figure 9) . Figure 10 – Contour plots of plated fracture callus; (A – C) fracture callus with a constant fracture gap and differ- ent conditions; (D – E) fracture callus with a constant loading condition and different fracture gaps; (E) wireframe representation of callus and plate positions. In Figure 11, posterior, medial and anterior cuts of the tibia at 150, 200, 250, 300 and 350 days show the stress shielding and damaged -based remodelling that occurs due to introdu ction of the plate. The remodelling framework considered cases with damage (D) and with no damage (ND). The highest degree of stress shielding occurred anteriorly ( Figure 11A) , while the least amount of stress shielding occurred posteriorly as seen in ( Figure 11C) . Stress concentrations were created locally, mainly in the interfacial regions between the titanium pins and cortical bone. Models containing damage (D) led to damage -based reabsorption around the pins, while the models without damage (ND) the hig h stresses resulted in positive SED remodelling around the pins. A significantly higher degree of stress shielding around the plate was observed for the models containing “no damage” than the models containing “damage” (Figure 11a) . 16 Figure 11 –Elastic modulus within cortical bone and fracture callus focus for plated tibia with 3mm fracture gap under LC(A) at 150, 200, 250, 300 and 350 days. Simulations were ran with damage (D) and with the no damage (ND). 4. Discussion In this study , a coupled biphasic model for bone fracture repair was developed that considers both fracture healing and remodelling phases of repair to investigate the performance of a fractured tibia in both non -plated and plated conditions. It was found that introduction of a titanium plate allowed the tibia to undergo successful healing at higher loading conditions and fracture gaps , compared to the non -plated versions. While these plated cases showed faster rates of repair in the healing phase, their performance was substantially different once they entered the remodelling phase, with all plated cases remodelling to a much lower ela stic modulus than any of the non -plated loading cases, with the coupled modelling framework predicting substantial regions of stress shielding. This 17 framework is one of the few implementations of coupled fracture healing and remodelling algorithms and incl uded several innovative approaches to smoothing, time -averaging and time incrementation in its implementation, thereby avoiding any unwanted abrupt changes between tissue phenotypes. This coupled model provides a better representation of tissue development in the fracture callus and cortical bone when compared with fracture healing alone and provides novel platform to investigate the long - term performance of orthopaedic fixation devices. This would include the more effective design of permanent fixators and provide an opportunity to optimise both spatial and temporal behaviour of bioabsorbable implants to limit unwanted effects of stress shielding. The coupled fracture repair model was shown to provide very good correlation to animal experiments carried out by Claes et al. (1995) , in terms of both predicted inter -fragmentary strain (see Figure 7) and the time sequence of the repair [54] . During the fracture healing phase, the simulations generally predicted the formation of the cartilaginous callus 2 -3 weeks from the start of the healing process [2]. The extremities of the fracture callus experience lower biphasic stimulus and stiffen more quickly . The bony bridge provides mechanical support lowering the biphasic stimuli at the fracture site, allowing for bone formation to occur in other regions such as near the bone marrow and fracture gap. The formation of a complete “bony bridge” was a pivotal aspect in achieving successful healing (see Figures 8-10), and the coupled fracture repair model used this feature as the criteria that determined whether the bone remodelling phase would be initiated (the time at which this occurred was not prescribed apriori ). Upon successful fracture healing, the non -plated models remodelled positively, tending towards the o riginal density of cortical bone (ρ = 1.649g/cm3, E = 17,000MPa). This was largely driven by the fact that a homeostatic load of 6.6MPa was assumed to be acting on the tibia after day 112 for all conditions, resulting in similar bone densities for the non -plated models. Due to the geometry of the non-plated models, they were remodelled in a symmetric manner, with the inner and outer callus undergoing negative remodelling, while the callus focus region underwent positive remodelling (see Figure 8). During al l transitions, the coupled fracture repair algorithm was implemented in such a way that abrupt changes in tissue modulus due to local stimuli were avoided. This was achieved through the implementation of several time incrementation controls and smoothing t echniques (see Section 2). For the latter, the change in elastic modulus due to the biphasic stimulus and SED were time averaged over 10 days. Therefore, a change in the elastic modulus of an element was due to its stimulus history, opposed to a single tim e instance as was commonly employed in other studies [55], [56], [57]. The implantation of a titanium plate played a critical role in healing and remodelling phases of bone repair, creating an uneven stress distribution across the fracture callus tha t resulted in substantially different healing outcomes between the plated and non -plated models. During fracture healing, regions directly adjacent to the plate experience lower biophysical stimuli than regions on the opposite side. Based on the mechano -regulatory framework, these elements transition through tissue phenotypes more quickly, leading to local stiffening, with some regions even being absorbed due to the lack of a biophysical stimulus during fracture healing (See Figure 10). Upon entering the re modelling phase, most of the stress is distributed through the implant, as it has a stiffness of approximately 6.5 times greater than cortical bone, thereby shielding the adjacent region. The proximity of the surrounding tissue to the titanium plate heavil y influenced the degree of stress -shielding. The highest degree of stress shielding occurred anteriorly, while the least amount of stress shielding occurred 18 posteriorly (see Figure 11). This subsequently resulted in the anterior portion of the “callus focu s” region to negatively remodel, while the posterior portion remodelled positively producing a cortical shell. Most of the bone loss around the titanium plate occurred between 0 – 6 months for the “damage” models with a significant reduction in bone loss f ollowing the 6 -months, matching experimental tests of an implanted titanium stem carried by Jaffe and Scott [58] . Interestingly, during the remodelling phases, the plated models also predicted further resorption in stress concentration regions where screws were fastened to the bone. The stress concentrations were generated due to (i) the dissimilar stiffness of titanium (E T = 110GPa) and cortical bone (E c = 17GPa) and (ii) the geometry of hole generated by the titanium pin. In this study, the role of dama ge-based reabsorption was during remodelling was considered . Within the “no -damage” models, the high stresses generated around the pins led to high levels of positive remodelling. This subsequently led to higher levels of stress being redistributed through the titanium plate increasing the stress -shielding effect. On the contrary, when damage was included in these models, these high stresses led to damage -based reabsorption in these regions (which is a completely separate effect to stress -shielding). Intere stingly, the effect of this damage -based resorption actually redistributed stress away from the pins and subsequently the plate, thereby substantially reducing the stress shielding effect. Damage -based reabsorption occurs much more rapidly than negative SE D remodelling, which may cause implant failure as the interfacial bone around the pins is reabsorbed, resulting in bone plate slippage. The healing performance for both the plated and non -plated simulation were directly impacted by the loading conditions (Figures 8 and 10) and fracture gap size (see Figure 9) . For the simulations in which the load was altered the fracture gap remained constant at 3mm, while for all different fracture gap sizes, the base loading condition LC(A) was used. Under different loa ding regimes (in the first 28 days), it was found that the time at which the “bony bridge” formed was different in both plated and non-plated scenarios. The plated models were able to successfully heal at much higher loads (e.g. LC(H) when compared with th e non -plated models (e.g. LC(C)). An increase in fracture gap size decreased the healing performance for both plated and non -plated simulations. Larger fracture gaps increased the biphasic stimulus induced within the fracture region for the same loading co nditions. For all fracture gaps the diffusion coefficient remained constant at 0.3mm/day. This, in turn, influenced the time taken for the callus focus region to become confluent in MSCs, and thus slowed the healing response. The simulated results obtained here mimic that of experimental results seen in a rat model conducted by Meeson et al (2019) , whereby a femoral osteotomy of different sizes (1mm, 1.5mm and 2mm) resulted in delayed healing over a 5 week period. The prediction of long -term outcomes of bone fracture repair represents an extremely challenging problem that involves highly complex interaction between physical, chemical and biological processes taking place at the fracture callus. This implies that the coupled healing and remodelling framework presented here is subject to several limitations. Both the fracture healing and remodelling algorithms rely on sets of phenomenological rules to drive cellular activity, according to a series of thresholds that dictate tissue formation at different stages of the process. For these, threshold parameters have been chosen based on existing experimental and computational literature across both algorithm types [26], [60] and this coupled framework uses a novel numeric al implementation, whereby abrupt changes between tissue phenotype is avoided through suitable smoothing, time -averaging and automatic time 19 incrementation strategies. In validating the coupled framework, the model predictions showed excellent agreement wit h inter -fragmentary movement (IFM) from an animal experiment, with differences below 0.1 mm seen throughout the healing phase. The model also produced typical patterns and appropriate time sequences that have been observed in both bone healing and bone rem odelling phases across animal and human studies [3], [54] . Therefore, this coupled model provides a suitable platform to investigate bone fracture repair and the long -term performance of orthopaedic fixation devices. 5. Conclusion In this study, a coupled biphasic model for bone fracture repair was developed that considers both fracture healing and remodelling phases of repair and investigates the performance of a fractured tibia in both non -plated and plated conditions. It was found that introduction of a titanium plate allowed the tibia to undergo successful healing at higher loading conditions and fracture gaps, compared to the non -plated versions, but resulted in substantial regions of stress shielding at later stages of repair. It was found healing tim e was increased when increases in either loading conditions (Figures 8 and 10) or fracture gap size (see Figure 9) were considered . This model representation of tissue development in the fracture callus and cortical bone when compared with fracture healing alone provides novel platform to investigate the long -term performance of orthopaedic fixation devices. This approach could be used to optimise the spatial and temporal load -bearing performance of bioabsorbable fixation devices, to enhance repair outcomes . Funding Funding support was provided by the Irish Research Council (IRC) Government of Ireland Postgraduate Scholarship (GOIPG/2017/ 2102). Conflicts of Interest/Competing interests The authors declare no conflict of interest(s) in the work presented. Code availability The model is implemented in the Abaqus/Standard finite element framework, with several user - defined field subroutines (USDFL D) written to predict fracture healing and bone remodelling. While these are not being made available through an open -source repository, the authors will consider requests for codes (or part thereof) through the corresponding author. 20