24 Numerical Simulation of Fracture Healing and Bone Remodeling
Besides biological factors, the local mechanical environment is known to play a crucial role in directing the fracture healing process. The remarkably strong dependency of osteogenesis on mechanical factors is likely to be the reason why engineering methods and computer models found their way into this specific field of orthopedic research rather early on. Within the last 20 years, different research groups have developed dynamic models in order to predict the time-dependent fracture healing processes. While early models mainly served to test the plausibility of existing tissue differentiation hypotheses, more complex and more thoroughly validated models have since been used more frequently to generate clinically relevant predictions and recommendations.
In any case, computer models are very helpful in advancing our understanding of complex systems, especially if such numerical models are tightly coupled with experimental research methods.
24.1 Understanding and Creating Computer Models of Biological Processes
Let us now address some basic simulation concepts one has to understand first, in order to be able to classify existing numerical models. Anyone building a model of any biological process has to answer the following fundamental questions first.
24.1.1 First Question: What Are the Main Players, the State Variables?
First, you should think about the main players that participate in your “game,” the so-called state variables or dependent variables. They describe the (changing) state of your system. They are the unknowns in your research question, the quantities you would like to measure, even if that is impossible experimentally. Start with one or only a few state variables and avoid unnecessary initial complexity; at first, include only the most important effects. It is your decision: You are the director; you can hire further players later on.
Examples of state variables:
Local osteoblast concentration in a fracture callus
Bone tissue distribution in a fracture callus
Bone mineral density distribution around a hip stem
24.1.2 Second Question: Do You Like to Get a Single Picture or a Movie? What Are the Independent Variables?
Next, you need to choose the independent variables (i.e., time and spatial coordinates). They awaken your players and provide them with room to live in; otherwise, they would just be constants.
Examples of independent variables:
Choose one, two, or three spatial coordinates (e.g., x, y, or z) in order to define the state variables as fields in a one-, two-, or three-dimensional space.
Add a time coordinate t if your problem is time-dependent and if you are interested in that time dependency. Introducing time means switching from a static to a dynamic approach.
By choosing the independent coordinates, you implicitly define the mathematical type of your problem and thereby the numerical method, the algorithm(s)—and thus, the kind of software—suitable for solving the model (see remarks in the next section).
24.1.3 Third Question: What Is the Story About? Name Effects and Processes
Now that you have your players listed, you are able to write down the interactions between them. Most dynamical models, especially those of biological processes, use the following quite general form of equations, which is much easier to read and to construct than you might think at first. Such a mathematical formalism likely fits your problem as well. For each of your previously chosen state variables, you need to devise an equation that describes how the state variable changes over time. This results in a system of equations of the following form:
Rate of change of variable 1 = Effect A + Effect B + …
Rate of change of variable 2 = Effect C + Effect D + …
This system describes how the state variables (left side) change over time depending on several effects (right side). The effects can increase or decrease the rate of change. They can depend on any of the state variables, even the one that can be found on the left side of the same equation, as well as on any of the independent variables, explicitly. Mathematically speaking, we are dealing with a (potentially coupled) set of (linear or nonlinear) first-order ordinary differential equations.
Example: Rate of change of bone mineral density = bone matrix production rate – bone absorption rate
The bone matrix production itself might be modeled as a function depending on the local osteoblast concentration. It would therefore be nice to have that osteoblast concentration as a further state variable available. To describe more complex effects, we need to also include derivatives of state variables with respect to the spatial coordinates, leading to partial differential equations.
24.2 Biological Processes of Fracture Healing
In order to be able to capture the underlying biological processes of fracture healing in a mathematical model, we first need to understand them in detail. Please refer to Chapter 9, “Biomechanics of Trabecular and Cortical Bone,” and Chapter 10, “Biomechanics of Fracture Fixation” for a detailed overview of the biological aspects of fracture healing and bone formation processes. In the context of this chapter, we only want to repeat briefly the most important aspects of the involved processes.
24.2.1 Processes of Fracture Healing
The natural form of fracture healing occurring in long bones is indirect or callus healing where the fragments, with or without artificial fixation, display a certain amount of so-called interfragmentary movement (IFM; i. e., relative movement toward each other).
The healing process (Fig. 24.1) starts with the phase of acute inflammation within the blood clot. All cells enclosed in the hematoma die within days after the trauma due to acute hypoxia. Leukocytes invade the area; granulocytes, histiocytes, and mastocytes start clearing the site of debris tissue. Then trauma and inflammatory reactions trigger the primary callus response: Underneath the periosteum, at some distance from the gap, intramembranous ossification starts and produces large amounts of woven bone. This primary response lasts for about 2 weeks.
In the meantime, revascularization of the hematoma begins and mesenchymal stem cells and fibroblasts invade the fracture site. The fibroblasts proliferate, produce collagen fibers, and replace the hematoma with well-vascularized granulation tissue. Within weeks, the initial granulation tissue is replaced by a cuff of fibrous connective tissue that encloses and connects both fracture fragments, forming a soft callus.
Near the fracture gap, where high strains predominate and prohibit both revascularization and intramembranous ossification, mesenchymal stem cells differentiate into chondroblasts that produce cartilage matrix. The produced cartilage further stabilizes the fracture site and provides the framework for the following ossification process.
Gradually, the fibrocartilage tissue mixture gets transformed into woven bone via endochondral ossification: Starting from the already existing woven bone that formed underneath the periosteum at some distance from the fracture gap, the ossification fronts of the two callus halves grow toward each other. While they approach, they increase their width as the mechanical strains in the vicinity of the fracture gap are far too high for the ossification process. This also increases the cross-sectional area and in turn the bending stiffness. Bony bridging through unification of the two callus wedges occurs typically in the periphery of the callus. When this happens, the fracture is considered to be “healed.” Besides the periosteal response, there is also an endosteal reaction, which is however less pronounced due to the delayed medullary “vascular response.”
24.2.2 Remodeling
Remodeling is a constantly ongoing process of resorption of old and formation of new bone to adapt the bone to new mechanical requirements, to repair microdamage, to regulate calcium homeostasis, to form the overall shape of the bone, or to transform woven into lamellar bone. Osteoclasts together with osteoblasts form a remodeling unit, where the preceding osteoclasts resorb bone via phagocytosis by secreting enzymes that break up the existing bone tissue. The following osteoblasts produce the extracellular matrix called osteoid that becomes mineralized, resulting in the formation of new bone tissue. This process is strongly controlled by the local mechanical stimuli.
24.2.3 Tissue Differentiation and Mechanoregulation Hypotheses
The discovery of the functional adaptation of living tissues to their daily mechanical loading goes back to the 19th century. Julius Wolff observed that the trabeculae in the femoral head apparently are aligned along the trajectories of principal mechanical stresses. He deduced in 1892 that local mechanical stimuli determine local tissue formation. 1 Wilhelm Roux further hypothesized that the development of specific tissues is associated with specific types of mechanical stresses. 2 Friedrich Pauwels improved this concept and hypothesized that formation of either cartilage or fibrous tissue would be stimulated by shape changing or volume changing strain states, respectively. 3 A further differentiation into bone tissue, however, requires a stable situation with only very low strains but no preference for a specific deformation type. Pauwels therefore concluded that the primary goal of a fracture treatment must be the stabilization of the fracture site such that bone tissue can form.
Beginning in the 1960s, Harold Frost explored the relations between mechanical loading and the remodeling process of bone tissue, but also other potentially mechanosensitive tissues over the course of more than four decades. 4 He suggested that the skeletal physiology—its local strength and architecture—is controlled by biological factors and mechanical influences in such a way that the mechanical strains always fall within a certain acceptable range (Utah paradigm of skeletal physiology). If the mechanical load exceeded an upper threshold, the tissue′s functional units would increase the strength of the tissue at that particular location. Respectively, if the stimulus dropped below a certain threshold, tissue would be resorbed until the local stimuli are again within the adapted window (Fig. 24.2).
This negative feedback loop (“mechanostat”) tends to plateau at equilibrium. Nearly all numerical remodeling simulations are based on this concept.
Jargon Simplified: Invariant
An invariant of the strain tensor (which can be represented by a matrix containing the strain components) is a value derived from the strain tensor′s components and that is independent from the choice of a particular coordinate system. Using strain invariants as mechanical stimuli is plausible, as the cells’ behavior should not depend on the arbitrary choice of a frame of reference.
Based on Pauwels ideas, Claes and Heigele 5 developed a first quantitative tissue differentiation function. Simon et al 6 enhanced that function and thermodynamically more consistently used only strain invariants as stimuli. Bone formation according to this function is possible only on existing bony surfaces with adequate blood supply and both strain stimuli lying within a certain range, called “medium” (Fig. 24.3).
Prendergast et al 7 used a more complex biphasic (poroelastic) material model to describe the biological tissues and introduced a tissue differentiation algorithm dependent on shear strain (distortional strain γ) in the solid phase and the fluid flow velocity in the interstitial fluid phase (Fig. 24.4). This function predicted bone formation for lower but not too low strains and fluid flows.
Though both hypotheses appear quite different at first sight, they both are able to explain the subsequent occurrence of fibrous soft tissue, cartilage, and bone, as both hypotheses associate tissue types with the magnitude of mechanical stimulation in an analogous order: Fibrous tissue can withstand high (in particular tensile) strains; chondrogenesis is the response to medium to high compressive and/or shear strains; and bone growth requires small but nonzero mechanical stimulation.