1,349
Views
0
CrossRef citations to date
0
Altmetric
Articles

Simulating localised cellular inflammation and substrate properties in a strain energy density based bone remodelling algorithm for use in modelling trauma

&
Pages 208-218 | Received 20 Sep 2017, Accepted 06 Feb 2018, Published online: 16 Feb 2018

Abstract

Bone responds to mechanical stimulus and a range of pre-existing finite element models have been suggested to reproduce the internal physiological structure of bone. Inflammation effects are not included in these models, yet inflammation is a key component of bone repair in trauma. Therefore, a model is proposed and tested here that extends these methods to include parameters that could be considered to represent the behaviour of bone remodelling when influenced by inflammation. The proposed model regulates remodelling based on findings from recent studies into the nature of heterotopic ossification, the formation of heterotopic bone, which have revealed information about the nature of bone after high levels of trauma. These parameters include consideration of the distance from the zone of trauma, the density of mesenchymal stem cells, and substrate stiffness as a trigger for cells becoming osteogenic. The method is tested on a two-dimensional plate model and shows that the new extended algorithm can produce a range of structures depending on inputs that could be used in the future to replicate physiological scenarios.

1. Introduction

Under normal conditions, the regulation of bone density is responsive to mechanical stimuli, with density increasing with higher bending stress. For example, this is found with resistance training (Layne and Nelson Citation1999). In addition, morphological changes are also found, for example, increased bone diameter tends to be present in the racquet forearm of tennis players (Ireland et al. Citation2013). Conversely, a reduction of bone density and strength is seen after bone underuse (Giangregorio and McCartney Citation2006) or stress shielding (Huiskes et al. Citation1987). Adaptive finite element modelling can function as a good approximation tool for showing the mechanical influence on bone remodelling (Mullender et al. Citation1994). These models are adaptive in the sense that the properties of the structure adapt or vary in response to feedback of results. The feedback can be considered as the error between some experienced stimulus and a steady state stimulus value. Previous authors have used a strain energy based modelling approach to model osteophytes (He and Xinghua Citation2006) and the development of ectopic bone formed after total disc replacement (Ganbat et al. Citation2014).

The behaviour of bone remodelling in response to inflammation may be dependent on other parameters. In order to understand how inflammation may affect bone remodelling, the clinical instance of heterotopic ossification (HO) acquired after tissue trauma can be examined. HO is the process of bone formation in tissues that are not usually osseous. The exact causes of HO are unknown, however, it is widely accepted that the inflammatory response of the body plays a major role. HO has been reported to occur in over 60% of injured service members returning from the recent conflicts in the Middle East (Alfieri et al. Citation2012) who suffered from traumatic and blast related injuries. Brown et al. Citation2010 reported the presence of HO in 20 out of 35 limbs in UK traumatic amputations and 134 out of 213 limbs in US amputations. Traumatically acquired HO is of current clinical interest and has such been the subject of a number of studies which have investigated the relationship between bone formation and trauma.

This paper aims to expand upon the strain energy based bone remodelling algorithm presented by Mullender et al. (Citation1994) to consider parameters that may reflect the remodelling of bone after injury. These parameters aim to reflect the effects of cellular inflammation whereby the likelihood of ossification is related to the localised level of inflammation, the level of strain in the tissue and the loading environment. In this study, the model will be compared against pre-existing tests on a 40 by 40 plane stress square plate in order to examine the effect of the introduced parameters through direct comparison with the existing model.

2. Mechanobiological theory

Within the square block, based on the existing theory presented by Mullender et al. (Citation1994), the chance of an element becoming stiffer is dependent on its loading and the loading in neighbouring elements. This theory can be extended to take into account cellular and inflammatory aspects of healing. It may be assumed that the results of inflammation may cause some cells to have a higher level of osteogenic precursor than other cells (or elements in this mathematical representation). After a traumatic injury, an inflammatory response will be triggered by the immune system. This involves increasing the blood flow in order to deliver necessary immune cells to the affected area, and to begin wound healing and regeneration processes. The migration of these cells and the effects of released cytokines and growth factors associated with inflammation increases the number of potentially osteogenic precursor cells in the area. This process may play a key role in bone modelling and remodelling. Furthermore, degrading bone fragment particles may initiate the proliferation and differentiation of osteoprogenitor cells (Nagi et al. Citation2002). The link between trauma and heterotopic bone formation is supported by a number of observations made in the case of traumatically acquired HO and in findings from studies on the behaviour of traumatically acquired HO (Potter et al. Citation2007; Jackson et al. Citation2009; Davis et al. Citation2011). Certain post wound treatments such as tissue debridement and irrigation have been noted to increase the probability of HO and this has been linked to the increased tissue trauma produced by these treatments (Potter et al. Citation2007). The application of positive and negative pressure dressings together with the use of antibiotics that aim to aid in wound healing may result in an over active tissue growth response (Ji et al. Citation2012) and tests have confirmed high proliferation of inflammatory and progenitor cells in patients at the site of future HO (Jackson et al. Citation2009; Davis et al. Citation2011). Davis et al. (Citation2011) hypothesised that the risk of HO increases with the number of progenitor cells whether they be osteogenic or not. Therefore, the level of an elements osteogenic precursor was considered to be related to the level of tissue trauma as well as the level of mechanical loading.

Potter et al. (Citation2007) found that amputations within the zone of injury are a risk factor for HO, therefore, the level of an elements osteogenic precursor was considered to increase with proximity to a defined site of trauma as well as with proximity to regions of higher mechanical loading. In the proposed algorithm, this can be reflected by a parameter that represents the level of mesenchymal stem cells within an element which is related to the location of the defined site of trauma or injury (Equation (Equation1)). In this study, only the effect of moving the point of trauma was examined to see how the algorithm affects the pattern of remodelling within elements.

The level of osteogenic precursor is not the only factor that can affect the chance of an element to stiffen. There must also be some consideration to the chance of osteogenic precursor actually resulting in osteogenesis. This may be highly dependent on how osteogenic or favourable the environment is. This refers to the physical environment and can be described by the natural substrate stiffness of tissue which is responsible for the transmission of forces to cells. Substrate stiffness affects stem cell lineage and can govern the chance of the progenitor cells differentiating into osteogenic cells (Engler et al. Citation2007; Zajac and Discher Citation2008; Cigognini et al. Citation2013). The assumption is made in this study that an element with a stiffness close to bone will have a higher probability of osteogenic differentiation and that those elements with a lower stiffness, such as that of soft tissue, will more likely result in myogenic differentiation (Clause et al. Citation2010). The differentiation of osteogenic cells is also governed by other factors such as systemic cytokines and hormones, pH levels and available oxygen. These parameters have been neglected from the model for simplicity in order to focus on the model to examine the mechanical and inflammatory influence on bone remodelling.

The effect of the parameters in the algorithm on the remodelling pattern was assessed by moving the location of high trauma or precursor to regions experiencing different levels of mechanical stimulus which thus create different osteogenic environments.

3. Mechanobiological algorithm of HO

A numeric representation of the level of osteogenic precursor is defined below as ρMSC, the density of mesenchymal stem cells.(1)

ρMSC is defined in Equation (Equation1) as a function of distance to injury point (y). The function f(y) is defined as a value between zero and 1, relative to the location of designated point of trauma or inflammation in the model space. For this study, the different functions for each scenario are shown in Table . ρMSC is limited between the maximum of 1 which represents the highest level of precursor and the minimum of 0.01. The minimum value should be small but greater than zero as even in healthy tissue there is a small amount of mesenchymal stem cells. This is supported by findings from a study by Davis et al. (Citation2011) where the number of progenitor cells was counted in uninjured tissue, injured tissues that did not develop HO, and injured tissues that did develop HO. The findings are tabulated in Table .

Table 1. Relative number of progenitor cells present in control, HO and non-HO sites based on data from Davis et al. (Citation2011).

In order to incorporate the chance of a cell becoming osteogenic due to substrate stiffness, the state of the osteogenic environment is described as a probability function based on element stiffness E.(2)

The function f(E) is defined in Table . The probability function P and the function of density of mesenchymal stem cells ρMSC can be combined to create a single factor that regulates remodelling of an element in a physiologically realistic way. This factor shall be named .

Table 2. Parameters used in the test of the novel model for HO growth.

Equation (Equation3) shows the suggested relationship among and β. α is a constant and is equal to 1.5(3)

The derivation of Equation (Equation3) was found as follows. In order to gain insight into how the location and level of inflammation may affect the likelihood of an element stiffening, the instance of HO was studied. In soft tissue, or tissue with a stiffness much lower than bone (and thus a low P), HO will not occur unless ρMSC is high. Even if ρMSC is high, the amount of remodelling should be lower than in normal bone remodelling regions as new bone can only form if the environment permits it. Therefore, at high levels of precursor (ρMSC), the overall chance of ossification (and the value of β) becomes more dependent on the level of substrate stiffness P (a function of stiffness E) when P is decreased to low levels. At low levels of precursor (ρMSC), β is more sensitive to changes in on stiffer substrates like bone callus and cartilage than soft tissue. This is because bone remodelling may still occur but at variable rates depending on the available precursor. On soft tissue however, low to mid-range levels of precursor is much less likely to result in bone modelling. In bone, regardless of the amount of inducing agent, normal bone remodelling should apply. In this case, the same algorithm as in Mullender et al. (Citation1994) applies and β should be close or equal to 1.

The proposed power law relationship (Equation (Equation3)) satisfies the condition that when substrate environment P is 1 (in the case of bone), β will also be 1, independent of the level of precursor ρMSC. It also satisfies the condition that β is more sensitive to changes in ρMSC when the level of substrate stiffness is high. The constant α of 1.5 in Equation (Equation3) was chosen in order to produce a general trend in β that reflected the physiological behaviour described above. Hence, the value was chosen by determining what value results in the sensitivity of β becoming greater to changes in P at low levels, whilst at high levels of ρMSC. High was considered as a value higher than the level of progenitors found in injured tissues that did not develop HO (>0.70 taken from Table ). Figure shows how the change in β increases when substrate stiffness levels decrease for a constant value ≥ 1.4. A value of 1.5 was selected as it offered a smoother trend than a value of 1.4 whilst still minimising the value of β at low substrate stiffness levels. Figure displays the suggested relationship among and the resulting β,  and highlights some real life scenarios and where the assigned values for these scenarios lie. Table summarises these scenarios giving a brief explanation for each of the suggested levels for the associated parameters and β.

Figure 1. Sensitivity of regulatory factor to level of substrate stiffness when the level of precursor is relatively high. Each curve shows the calculated change in in relation to the change in for different values of the constant α (1 to 1.6).

Figure 1. Sensitivity of regulatory factor to level of substrate stiffness when the level of precursor is relatively high. Each curve shows the calculated change in in relation to the change in for different values of the constant α (1 to 1.6).

Figure 2. Relationship between and . Dashed rectangle represents region of normal bone healing, Oval represents region of callus healing, shaded rectangle represents that of aberrant wound healing, dotted rectangle encompasses values for soft tissue healing and pentagon covers region of activity in stable soft tissue.

Figure 2. Relationship between and . Dashed rectangle represents region of normal bone healing, Oval represents region of callus healing, shaded rectangle represents that of aberrant wound healing, dotted rectangle encompasses values for soft tissue healing and pentagon covers region of activity in stable soft tissue.

β, can be added to Equation (6) in the paper by Mullender et al. (Citation1994) to produce the following equation and regulate the change in density:(4)

In Equation (Equation4), τ is a time constant and ϕ is the level of strain energy based stimulus defined in Equation (Equation5) (Mullender et al. Citation1994). τ is kept at a value of 0.1. This value ensures that the change in density is just over 20% of the total range per iteration. This results in a visually steady rate of change as opposed to large jumps in the rate of change of density.(5)

Equation (Equation5) describes that the stimulus at position x and time t is the sum of the error between the experienced strain energy and the steady state reference stimulus K for every element (from 1 to N) where N is 1600, the total number of elements. In this model, each element is considered to behave as a region of tissue with an individual sensor point in the centre. The exponent is a spatial influence parameter considering a sensor influence factor D and the distance between each sensor to position x. In this study, the reference stimulus K was taken to be 0.25 J/g as per the study by Mullender et al. (Citation1994). The sensor influence factor D was also kept to be the same as per the Mullender et al. (Citation1994) study which is the distance between the central nodes of adjacent elements.

4. Testing the novel algorithm to predict shape of HO in a standard test specimen – a 2D plate

4.1. Aim

The model presented in Mullender et al. (Citation1994) reproduced a load dependant discontinuous trabecular architecture in a 1 mm2 square plate. The aim of this study is to apply the new remodelling algorithm (Equation (Equation4)) to the 2D square plate used in the paper by Mullender et al. (Citation1994). This allows for a standardised test of the newly developed model which will highlight how it behaves differently. The simplicity of the loaded plate mesh enables the effects of the loading on the model to be clearly visualised.

4.2. Method

The function of β is to maximise bone remodelling in areas likely to behave as bone, and minimise it in regions that should not change their state.

The value of ρMSC (the quantity, or density, of mesenchymal stem cells) reduces with distance from the injury point (Equation (Equation1)). Six different hypothetical injury scenarios were simulated (labelled Scenarios 1–6). These are shown in Figure (a–e). The patterns shown are the resulting density distribution from that which is obtained when using the remodelling algorithm presented in Mullender et al. (Citation1994) publication with the same ramped compressive load on a 40 by 40 plate (Figure ).The white spot in these figures shows where ρMSC is at its maximum of 1, hence representing a high risk region for element stiffening (such as in the case of high levels of precursor increasing the risk of HO). ρMSC reduces linearly in the radial direction reaching a minimum at the displayed outer rings. For Scenarios 1–5, the minimum value of ρMSC is 0.01 (Table 2). For Scenario 6 the outer ring marks a minimum ρMSC of 0.5. Outside of this ρMSC it immediately drops to 0.01. The different scenarios were chosen examine how the different parts of the proposed algorithm alter the remodelling. The effect of changing the function ρMSC is assessed by the comparison between Scenarios 1–3 and Scenarios 4–6. Scenarios 1–3 have a gradually varying value of ρMSC throughout the area of the model. Scenarios 4–6 have localised areas of high precursor ρMSC (shown by the circles overlaid on Figure and described by the equations in Table ), outside of which the value of ρMSC is negligible. The effect of the function P is assessed by moving the location of the high risk HO sites relative to the surrounding environment. Scenario 1 tested the remodelling when the peak value of mesenchymal stem cells was located in the middle of a high stiffness strut, representing an area of naturally high substrate stiffness such as bone. Scenario 2 tested the remodelling when the peak value of mesenchymal stem cells was located in an area of low stiffness surrounded by high stiffness struts, such as in a soft tissue region adjacent to bone. Scenario 3 tested the remodelling when the peak value of mesenchymal stem cells was located in a region with little stimulus and low stiffness.

Figure 3. Density results recreated from the method used in the paper by Mullender et al. (Citation1994) after applying a ramp load increasing in magnitude from right to left.

Figure 3. Density results recreated from the method used in the paper by Mullender et al. (Citation1994) after applying a ramp load increasing in magnitude from right to left.

Table 3. Healing scenarios and estimated values for regulating parameters and β.

Figure 4. Depiction of simulated Scenarios 1–6. The minimum value of the function of is overlaid as a ring on the density distribution found when no cellular effects are employed ( = 1). For Scenarios 1–5, the minimum value of is 0.01 (Table ). For Scenario 6 the outer ring marks a minimum of 0.5. Outside of this it immediately drops to 0.01. The white dot in the centre of the rings shows where is at its maximum of 1. Hence, reduces radially from the central white dot to the black circular outline. The different scenarios analyse the effects of the proposed algorithm by moving the hypothetical level of precursor () to different areas with respect to the surrounding stiffness levels of the elements and by altering the function of itself.

Figure 4. Depiction of simulated Scenarios 1–6. The minimum value of the function of is overlaid as a ring on the density distribution found when no cellular effects are employed ( = 1). For Scenarios 1–5, the minimum value of is 0.01 (Table 2). For Scenario 6 the outer ring marks a minimum of 0.5. Outside of this it immediately drops to 0.01. The white dot in the centre of the rings shows where is at its maximum of 1. Hence, reduces radially from the central white dot to the black circular outline. The different scenarios analyse the effects of the proposed algorithm by moving the hypothetical level of precursor () to different areas with respect to the surrounding stiffness levels of the elements and by altering the function of itself.

The regions of the plate in Figure (d) and (e) outside of the ring marking the minimum of ρMSC represent control tissue (Table ). The low value of ρMSC prevents any remodelling or change in density. Hence the material does not change its state. Aberrant wound healing or callus environment is represented by the regions in the scenarios where ρMSC is high.

The equations that describe the function of ρMSC for each scenario at position (xy) where x and y are the coordinates of the sensor located in the element are shown in Table .

The probability function P was based on the final density distribution of the 40 by 40 plate model in Mullender et al. (Citation1994). The maximum Young’s modulus was 303 MPa. The modulus of the greyer areas was in the range of 100 to 300 MPa. The probability function was set as follows:

The above is described in equation form in Table for the parameter P along with a complete summary of the parameters used for each case. The initial density for every element was set as 0.8 g/cm3 which follows from the Mullender et al. (Citation1994) publication. The plate was loaded with a ramped pressure decreasing from 10 N/mm2 at the left to 0 N/mm2 on the right.

There is no finite end point identified in the simulations. The end point of the simulation was chosen to be at the point where the change in mass reached a plateau and no visual changes in the structure could be observed. A plateau was consistently seen by iteration 8000 and thus results will be shown at iteration 8000.

5. Results

This section displays the findings of each test. Figure shows the density distribution at iteration 8000 for Scenarios 1–6 demonstrating the sensitivity of the final geometry to location and spread of the modelled trauma. The convergence of the mass of the system for Scenario 3 is shown in Figure , demonstrating mass convergence by iteration 8000. Figure shows the initial values for Scenarios 1–3 and the initial level of stimulus . Figure shows the progression of distribution for Scenario 3 results demonstrating initial general densification, followed by prime strut formation and then smaller strut formation and Figure plots the mean absolute difference between and for Scenario 3 over the course of the simulation. Figure shows the density distribution at iteration 8000 after cellular effects are removed from iteration 3680 which is when the mean absolute difference between and drops below 0.02 (Figure ). This demonstrates that cellular effect are dominant initially and are then negligible after the structure has started to form.

Figure 5. Density distribution from the novel algorithm for HO growth Scenarios 1–6. (a)–(c) show the response when the peak level of precursor is moved to regions experiencing different levels of stimulus. The resulting structure is similar to those seen under normal bone remodelling, however, an extra number of seemingly less functional elements appear on the right. (d)–(f) shows the response when the level of precursor to HO reaches a minimum within the remodelling plate area. The result is confined remodelling. The difference between (d) and (f) is a higher level of precursor within the remodelling elements for (f).

Figure 5. Density distribution from the novel algorithm for HO growth Scenarios 1–6. (a)–(c) show the response when the peak level of precursor is moved to regions experiencing different levels of stimulus. The resulting structure is similar to those seen under normal bone remodelling, however, an extra number of seemingly less functional elements appear on the right. (d)–(f) shows the response when the level of precursor to HO reaches a minimum within the remodelling plate area. The result is confined remodelling. The difference between (d) and (f) is a higher level of precursor within the remodelling elements for (f).

Figure 6. Convergence of mass of Scenario 3.

Figure 6. Convergence of mass of Scenario 3.

Figure 7. Mechanical stimulus (top) and regulatory factor (bottom) values for the first iteration of Scenarios 1, 2 and 3 of the novel algorithm for HO growth. Peak is considerably further from peak in Scenario 3 in comparison to 1 and 2.

Figure 7. Mechanical stimulus (top) and regulatory factor (bottom) values for the first iteration of Scenarios 1, 2 and 3 of the novel algorithm for HO growth. Peak is considerably further from peak in Scenario 3 in comparison to 1 and 2.

Figure 8. Progression of in Scenario 3 of the novel algorithm for HO growth. In early iterations the distribution of is dominated by the distribution of the level of precursor. In later iterations the distribution of is dominated by the level of mechanical stimulus.

Figure 8. Progression of in Scenario 3 of the novel algorithm for HO growth. In early iterations the distribution of is dominated by the distribution of the level of precursor. In later iterations the distribution of is dominated by the level of mechanical stimulus.

Figure 9. Difference between and for Scenario 3 over time.

Figure 9. Difference between and for Scenario 3 over time.

Figure 10. Density distribution for Scenario 3 after cellular effects are disabled from iteration 3680. The results are similar to when cellular effects are included throughout because of their diminishing effect over time.

Figure 10. Density distribution for Scenario 3 after cellular effects are disabled from iteration 3680. The results are similar to when cellular effects are included throughout because of their diminishing effect over time.

6. Discussion

The results show that the newly developed method can produce results similar to the previous methods as shown in Mullender et al. (Citation1994). However, the different scenarios resulted in different structures. The effect of the function can be seen to allow for control of remodelling regions (Figure (d)–(f)).

The effect of locating the trauma point over areas experiencing different amounts of stimulus was also seen to cause different remodelling responses (Figure (a)–(c)). Adjusting the function of ρMSC was also seen to locally affect remodelling (Figure (f)). Scenarios 1 and 2 produced results which most closely matched those in the literature (Mullender et al. Citation1994). However, extra narrow columns of bone appear on the right hand side of the plate where loading is smaller. The narrow struts in the original paper by Mullender et al. branch to the denser thicker struts which are located near the fixation point of the model. The columns appearing in the current study do not connect to a denser supporting strut. These struts do not seem to be functionally supporting load and may indicate a slightly less optimised and disorganised structure. This is reasonable due to the part of the algorithm that increases element stiffening due to hypothetical inflammation rather than loading.

Scenario 3, however, resulted in a different structure consisting of one thicker branch of stiffened material at the base. In the case of Scenario 3, and consequently β is increased in areas which before had a low value of remodelling stimulus ϕ. This seems to result in a different density distribution and remodelling path in comparison to when β is set to be high in regions where ϕ is already high as in the case of Scenarios 1 and 2. Figure shows the initial β values for Scenarios 1, 2 and 3 and the initial ϕ results. This suggests that when increased remodelling occurs in an unusual area the structure may be more significantly altered.

Figure depicts the progression of β distribution for Scenario 3. The function of β is initially dominated by the trauma location but over time becomes dominated by the current structure. Hence, initially greater remodelling or activity will occur in a high trauma area but later, effects of remodelling apply to the denser bonelike areas whereas reduced effects of remodelling apply to the less dense non-bone areas. At a certain point in time, the effects of become insignificant. This is supported by Figure which demonstrates the decreasing difference between β and the probability function P which is defined by element stiffness. The close similarity between Figures and c further depicts that once β becomes dominated by the structure, removing the effects of cellular activity has little effect on the progression of remodelling. This can be clinically related to the fact that inflammation and cellular effects from healing processes diminish over time, however, bone remodelling is continuous. The initial formation of HO may be heavily associated with inflammation but the continual and more gradual change in bone geometry is due to the nature of bone remodelling in response to local loading. In a study by Nagi et al. (Citation2002), HO was triggered in 6 out of 7 patients who experienced a shoulder dislocation within the first 7 days of shoulder arthroplasty surgery but not when dislocation occurred later. This suggests that effects from trauma are only influential whilst the inflammatory response is active. This is reiterated by the fact that prophylaxis of HO such as non-steroidal anti-inflammatory drugs and radiation therapy that function by inhibiting inflammatory responsive mechanisms and newly proliferating and differentiation osteoprogenitor cells, need to be administered early in treatment in order to be affective (Nauth et al. Citation2012). The effect of cyclic stretch increases mesenchymal stem cell proliferation in the first few days after in which the effect diminishes (Choi et al. Citation2007).

Scenarios 4–6 demonstrated that remodelling was confined to regions where was above the minimum. Outside of this, no density or material change occurs. Scenario 6 differed from Scenario 4 only in that was raised. The consequence of this was extended stiffening of elements in the locations closest to high loading. This is representative of a callus environment as stiffening occurs only in regions receiving mechanical stimulus within a localized area of high (Isaksson et al. Citation2009).

The ectopic bone models previously developed (He and Xinghua Citation2006; Ganbat et al. Citation2014) are appropriate for ectopic bone cases where mechanical stimuli are the main influence. It is accepted that osteophytes may form in response to the modified loading caused by degenerative changes in the vertebral discs. It has also been suggested in a clinical imaging study (Jin et al. Citation2013) that the morphology of ectopic bone after cervical disc replacement is dependent on the loading environment. However, in the case of traumatic injury, the complex morphology of heterotopic bone may be attributed to other factors. The current study produces different bone density distributions by incorporating additional parameters that intend to reflect the behaviour of cells which have been subjected to trauma and inflammation.

The main limitations of this study are related to the simplifications proposed. Validation is also difficult as the representative geometry is not physiological. Another simplification of the model is in that the model assumes a linear relationship between the rate of change of density and the stimulus value. The only limiting factor is the maximum and minimum density limits. However in reality, levels of stimulus and strain above a certain threshold are unfavourable for bone formation and can lead to resorption (Li et al. Citation2007). If the rate of repair is slower than the rate of micro damage, stress fractures can occur (García-Aznar et al. Citation2005). Further work could apply this remodelling algorithm to a physiologically representative finite element model such as a trans-femoral residual limb model that could be compared against radiographs of HO.

Supplementary material

Supplementary material relating to this article is available on request. Please contact [email protected] or [email protected].

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

This work was supported by the Engineering and Physical Sciences Research Council, [grant number EP/L504786/1] and was delivered with the Royal British Legion Centre for Blast Injury Studies.

Acknowledgements

The authors would like to thank Andy Bell, Senior Technical Consultant at MSC Software, Hampton, Greater London, United Kingdom for technical support in the Marc Mentat Software.

References

  • Alfieri KA, Forsberg JA, Potter BK. 2012. Blast injuries and heterotopic ossification. Bone & Joint Res. 1(8):192–197.
  • Brown KV, Dharm-Datta S, Potter BK, Etherington J, Mistlin A, Hsu JR, Clasper JC. 2010. Comparison of development of heterotopic ossification in injured US and UK armed services personnel with combat-related amputations: preliminary findings and hypotheses regarding causality. J Trauma. 69(Supplement):S116–S122.10.1097/TA.0b013e3181e44cc7
  • Choi K-M, Seo Y-K, Yoon H-H, Song K-Y, Kwon S-Y, Lee H-S, Park J-K. 2007. Effects of mechanical stimulation on the proliferation of bone marrow-derived human mesenchymal stem cells. Biotechnol Bioprocess Eng. 12(6):601–609.10.1007/BF02931075
  • Cigognini D, Lomas A, Kumar P, Satyam A, English A, Azeem A, Pandit A, Zeugolis D. 2013. Engineering in vitro microenvironments for cell based therapies and drug discovery. Drug Discovery Today. 18(21–22):1099–1108.10.1016/j.drudis.2013.06.007
  • Clause KC, Liu LJ, Tobita K. 2010. Directed stem cell differentiation: the role of physical forces. Cell CommAdhes. 17(2):48–54.10.3109/15419061.2010.492535
  • Davis TA, OʼBrien FP, Anam K, Grijalva S, Potter BK, Elster EA. 2011. Heterotopic ossification in complex orthopaedic combat wounds: quantification and characterization of osteogenic precursor cell activity in traumatized muscle. J Bone Joint Surg Am. 93(12):1122–1131.10.2106/JBJS.J.01417
  • Engler AJ, Rehfeldt F, Sen S, Discher DE. 2007. Microtissue elasticity: measurements by atomic force microscopy and its influence on cell differentiation. Methods Cell Biol. 83:521–545.10.1016/S0091-679X(07)83022-6
  • Ganbat D, Kim K, Jin YJ, Kim YH. 2014. Heterotopic ossification in cervical total disk replacement: a finite element analysis. Proc Inst Mech Eng H. 228(2):200–205.10.1177/0954411914522024
  • García-Aznar JM, Rueberg T, Doblare M. 2005. A bone remodelling model coupling microdamage growth and repair by 3D BMU-activity. Biomech Model Mechanobiol. 4(2-3):147–167.10.1007/s10237-005-0067-x
  • Giangregorio L, McCartney N. 2006. Bone loss and muscle atrophy in spinal cord injury: epidemiology, fracture prediction, and rehabilitation strategies. J Spinal Cord Med. 29(5):489–500.10.1080/10790268.2006.11753898
  • He G, Xinghua Z. 2006. The numerical simulation of osteophyte formation on the edge of the vertebral body using quantitative bone remodeling theory. Joint Bone Spine. 73(1):95–101.10.1016/j.jbspin.2005.03.019
  • Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Slooff TJ. 1987. Adaptive bone-remodeling theory applied to prosthetic-design analysis. J Biomech. 20(11-12):1135–1150.10.1016/0021-9290(87)90030-3
  • Ireland A, Maden-Wilkinson T, McPhee J, Cooke K, Narici M, Degens H, Rittweger J. 2013. Upper limb muscle–Bone Asymmetries and Bone Adaptation in Elite Youth Tennis Players. Med Sci Sports Exercise. 45(9):1749–1758.10.1249/MSS.0b013e31828f882f
  • Isaksson H, Gröngröft I, Wilson W, van Donkelaar CC, van Rietbergen B, Tami A, Huiskes R, Ito K. 2009. Remodeling of fracture callus in mice is consistent with mechanical loading and bone remodeling theory. J Orthop Res. 27(5):664–672.10.1002/jor.v27:5
  • Jackson WM, Aragon AB, Bulken-Hoover JD, Nesti LJ, Tuan RS. 2009. Putative heterotopic ossification progenitor cells derived from traumatized muscle. J Orthop Res. 27(12):1645–1651.10.1002/jor.20924
  • Ji Y, Christopherson GT, Kluk MW, Amrani O, Jackson WM, Nesti LJ. 2012. Heterotopic ossification following musculoskeletal trauma: modeling stem and progenitor cells in their microenvironment. Adv Exp Med Biol. 720:39–50.10.1007/978-1-4614-0254-1
  • Jin YJ, Park SB, Kim MJ, Kim KJ, Kim HJ. 2013. An analysis of heterotopic ossification in cervical disc arthroplasty: a novel morphologic classification of an ossified mass. Spine J. 13(4):408–420.10.1016/j.spinee.2012.11.048
  • Layne JE, Nelson ME. 1999. The effects of progressive resistance training on bone density: a review. Med Sci Sports Exercise. 31(1):25–30.10.1097/00005768-199901000-00006
  • Li J, Li H, Shi L, Fok AS, Ucer C, Devlin H, Horner K, Silikas N. 2007. A mathematical model for simulating the bone remodeling process under mechanical stimulus. Dent Mater. 23(9):1073–1078.10.1016/j.dental.2006.10.004
  • Mullender MG, Huiskes R, Weinans H. 1994. A physiological approach to the simulation of bone remodeling as a self-organizational control process. J Biomech. 27(11):1389–1394.10.1016/0021-9290(94)90049-3
  • Nagi O, Dhillon M, Batth H. 2002. Heterotopic ossification after total hip arthroplasty: a review of etiopathogenesis, risk factors and treatment modalities. Ind J Orthop. 36(4):225–233.
  • Nauth A, Giles E, Potter BK, Nesti LJ, OʼBrien FP, Bosse MJ, Anglen JO, Mehta S, Ahn J, Miclau T, Schemitsch EH. 2012. Heterotopic ossification in orthopaedic trauma. J Orthop Trauma. 26(12):684–688.10.1097/BOT.0b013e3182724624
  • Potter BK, Burns TC, Lacap AP, Granville RR, Gajewski DA. 2007. Heterotopic ossification following traumatic and combat-related amputations. Prevalence, risk factors, and preliminary results of excision. J Bone Joint Surg. 89(3):476–486.
  • Zajac A, Discher DE. 2008. Cell differentiation through tissue elasticity-coupled, myosin-driven remodeling. Curr Opin Cell Biol. 20(6):609–615.10.1016/j.ceb.2008.09.006