http://orcid.org/0000-0002-2278-2934
Figures
Abstract
Spinal fusion is a standard surgical treatment for patients suffering from low back pain attributed to disc degeneration. However, results are somewhat variable and unpredictable. With fusion the kinematic behaviour of the spine is altered. Fusion and/or stabilizing implants carrying considerable load and prevent rotation of the fused segments. Associated with these changes, a risk for accelerated disc degeneration at the adjacent levels to fusion has been demonstrated. However, there is yet no method to predict the effect of fusion surgery on the adjacent tissue levels, i.e. bone and disc. The aim of this study was to develop a coupled and patient-specific mechanoregulated model to predict disc generation and changes in bone density after spinal fusion and to validate the results relative to patient follow-up data. To do so, a multiscale disc mechanoregulation adaptation framework was developed and coupled with a previously developed bone remodelling algorithm. This made it possible to determine extra cellular matrix changes in the intervertebral disc and bone density changes simultaneously based on changes in loading due to fusion surgery. It was shown that for 10 cases the predicted change in bone density and degeneration grade conforms reasonable well to clinical follow-up data. This approach helps us to understand the effect of surgical intervention on the adjacent tissue remodelling. Thereby, providing the first insight for a spine surgeon as to which patient could potentially be treated successfully by spinal fusion and in which patient has a high risk for adjacent tissue changes.
Citation: Rijsbergen Mv, van Rietbergen B, Barthelemy V, Eltes P, Lazáry Á, Lacroix D, et al. (2018) Comparison of patient-specific computational models vs. clinical follow-up, for adjacent segment disc degeneration and bone remodelling after spinal fusion. PLoS ONE 13(8): e0200899. https://doi.org/10.1371/journal.pone.0200899
Editor: María Angeles Pérez, Universidad de Zaragoza, SPAIN
Received: October 30, 2017; Accepted: July 4, 2018; Published: August 30, 2018
Copyright: © 2018 Rijsbergen et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The data is owned by researchers who are not authors on the current study and belong to groups/institutions that formed the EU FP7 MySpine consortium. Data can be made available to qualified researchers upon request to the data steward, Dragan Bošnački (Email: d.bosnacki@tue.nl; tel: +31402475159). The data is kept at the different institutes according to project regulations and the institutions involved follow national regulations with regard to data storage. The installation of the data storage policy is under development. Until this is finished, data will be securely backed up on local disk clusters for at least 10 years.
Funding: The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7-ICT-2009-6) under grant agreement n° 269909, http://cordis.europa.eu/project/rcn/97394_en.html.
Competing interests: The authors have declared that no competing interests exist.
Abbreviations: AF, Annulus Fibrosis; BC, Boundary Condition; BV/TV, Bone density (volume fraction); BMD, Bone Mineral Density; CS, Chondroitin Sulfate; ECM, Extra Cellular Matrix; FCD, Fixed Charge Density; FU, Follow-Up; FE, Finite element; GAG, GlycosAminGlycans; IP, Integration Point; IVD, InterVertebral Disc; KS, Keratan Sulfate; MDC, Mean Detectable Change; MS, mechanical stimulus; MW, Molecular Weight; NP, Nucleus Pulposus; SED, Strain Energy Density
Introduction
Almost everyone experiences low back pain during their lifetime [1,2]. This high prevalence and its relative morbidity contribute to its large socio-economic burden [2]. Although the exact cause of low back pain is often unclear, there is a strong association with intervertebral disc degeneration [3,4]. If conservative symptomatic treatment fails, surgical intervention is considered [5]. Spinal fusion is the standard surgical treatment, however, results are somewhat variable and unpredictable [6,7]. Moreover, it has been demonstrated that fusion is associated with a risk for accelerated disc degeneration at neighbouring levels, most likely due to the loss in flexibility at the fused level causing altered loading conditions at neighbouring levels [8,9,10,11]. Presently, there is some ambiguity in the literature concerning whether fusion is the cause for accelerated disc degeneration at adjacent levels. Some studies report accelerated disc degeneration while others do not (see e.g. [11]). A possible explanation for these conflicting conclusions is that the outcome of each study is strongly affected by patient-specific factors, such as spinal geometry, level of activity, osteoporosis, etc. This, in turn, would suggest that a patient-specific evaluation is needed to better select fusion patients.
Whereas fusion may directly change the loading conditions at neighbouring discs, it may also change loading conditions at neighbouring vertebrae. According to Wolff’s law, changes in loading will lead to changes in bone density and architecture, which, in turn, may also lead to local changes in disc loading at neighbouring levels. However, similar to the disc, also the reported findings of changes in bone mineral density (BMD) do not conform. Some researchers reported a decrease in BMD at the adjacent level [12] while others found an increase [13,14]. Possibly these differences reflect differences in loading of the vertebrae before and after the surgical intervention. This would suggest that changes in loading may be non-uniform, and, again, that a patient-specific evaluation is needed to better predict such changes.
The most advanced type of patient-specific models available nowadays are computer models based on clinical images. For the disc, several such models have been presented in the literature [15,16,17,18,19]. These models, however, typically represent only a specific stage of degeneration [20,15,16,17,18,19,21], mostly in a phenomenological manner and cannot account for continuous changes in the biochemical composition and structure of the disc over time during degeneration. Recently, a composition based FE model of an Inter Vertebral Disc (IVD) was developed [22]. In this model, the behaviour of the IVD is modelled by relating the local extra cellular matrix (ECM) composition and organization directly to its mechanical behaviour. Thus, the annulus fibrosus (AF) and nucleus pulposus (NP) properties depend directly on its constituent (water, fixed charged density, collagen and ground substance) material properties proportional to their content within the tissues. However, in order to become a predictive model, it should also account for changes in biochemical composition over time due to changes in cell activity. Several models have been developed to predict such changes in connective tissue cell activity and ECM based on mechanical loading conditions sensed by the cells [23]. Although such models have been mainly used and validated to explain fracture healing [24], they in general can also predict tissue differentiation processes in other situations. Presently, however, no model is available that couples such advanced composition-dependent constitutive models of the disc with mechanoregulated tissue differentiation models.
For bone, several remodelling theories based on Wolff’s law have been formulated to predict changes in bone density related to changes in mechanical loading [25]. It has been demonstrated that such models can successfully predict changes in bone density after total hip [26] and knee arthroplasty [27]. In the case of spinal fusion, however, the situation is more complicated since changes in overall loading for neighbouring segments are expected due to the reduced flexibility of the fused segment, but also changes in local load transfer due to degeneration of adjacent discs. To study such changes, a coupled analysis of bone remodelling and disc degeneration is required in combination with a whole spine model to predict changes in overall loading.
The aim of this study was to develop such a coupled and patient-specific mechanoregulated model to predict disc generation and changes in bone density after spinal fusion and to validate the results relative to patient follow-up data. Therefore, first, the constituent based FE model of an IVD, including vertebrae and ligaments [22], was extended to a mechanoregulated IVD model and it was demonstrated that a stable (steady-state) healthy IVD with proper tissue morphology was maintained under physiological loading. Next, the IVD model was extended with a previously developed bone remodelling theory [28] and extended to a patient-specific (pt-specific) algorithm. This was done by coupling it to a full lumbar spine FE model [29] as was developed in the framework of the EU-funded project MySpine. From this full lumbar spine model, pt-specific changes in loading due to fusion at the adjacent level were derived and subsequently, converted to the IVD model. By doing this, tissue adaptation at the adjacent level was simulated for both bone and IVD tissue. To demonstrate the clinical feasibility of this approach, patient-specific models were made for 10 patients that underwent fusion surgery and the predicted computational results were compared to clinical follow-up data.
Material and method
First, the constituent based FE model of an IVD was extended to a mechanoregulated IVD model. To do so, a mechanoregulated tissue differentiation theory [23] was implemented, resulting in a stable healthy IVD that serves as starting point for IVD tissue adaptation simulations.
Disc adaptation framework
Steady state tissue composition.
The constituent based FE model of a healthy IVD [22] with generalized geometry was used as starting point (Fig 1). In this model, the NP and AF tissue properties are a function of their biochemical composition and microstructural organization and the mechanical behaviour is described according to [22] (Eq 1):
(Eq 1)
where ns,0 is the initial solid volume fraction, J the determinant of the deformation tensor F,σnf the stress in the non-fibrillar matrix,
the isotropic stress in the collagen fibres,
the tensile stress in the ith fibril,
the fibril density, totf the amount of fibrils, μf the water chemical potential, I the unit tensor and Δπ the osmotic pressure relative to the external physiological salt concentration (calculated based on the fixed charge density FCD). The distinction between IVD tissue properties, i.e. NP vs. AF, results from the differences in the distribution and organisation of the extracellular constituents (Table 1). For more detail, see [22].
To the generic IVD model a general daily loading pattern was applied, mimicking a loading cycle of a healthy generic person. Based on this loading pattern, per integration point (IP) a deviatoric shear strain and fluid velocity was obtained. These serve as input for the mechanoregulation algorithm and based on this algorithm, per IP, a preferred tissue phenotype is determined. The ECM content is adapted to this preferred phenotype with a fixed time step. After this adaptation step, the deviatoric shear strain and fluid velocity are re-determined and the whole cycle is repeated until no changes in ECM content were simulated.
Biochemical input parameters, for the IVD model [22]. The biochemical parameters are only used during the initial simulation of the disc adaptation. Thereafter, the mechanical stimulus determines per integration point what the content is.
To simulate IVD tissue remodeling, the model was extended by defining the cells and corresponding ECM components per integration point (IP). Two types of cells were assumed to be present: fibroblasts and chondrocytes. The fibroblasts could initially be mainly found within the AF tissue, while the chondrocytes are responsible for maintaining the NP. Both types of cells maintain their respective phenotypic ECM, consisting of collagen and glycosaminoglycans (GAG). Collagen content is a direct input for the disc model [22], whereas GAG content is converted into fixed charge density (FCD). This was done by taking into account the molecular weight (MW) and negative charges per molecule of chondroitin sulfate (CS) and keratan sulfate (KS) present in fixed proportions as GAG molecules [30], as well as the fluid fraction (volume based) at each IP. Based on this new matrix content (collagen and GAG), the model was allowed to come to swelling equilibrium and thus, allowed to come to its preferred water content (self-determined).
The behaviour of each cell was determined by a mechanical stimulus (MS). The MS was determined by applying generic boundary conditions representing a normal daily loading cycle (axial compression, flexion, dynamic axial rotation and dynamic lateral bending, representing standing, sitting, periodic movements and walking respectively). Axial compression represents the action of external loads as well as muscle activation and is hence constantly present in almost all daily recreational and occupational activities. The magnitude of the compressive load (300 N and 30 N for day and night load, respectively) was based on in vivo studies during daily activities with upright posture [31,32,33,34,35]. In addition to this axial compressive load, three additional bending rotations [31,32,33,34,35] were applied (while the axial compressive load was still present): flexion (2°) in combination with a creep period of half an hour to represent sitting; dynamic lateral bending (1.5°) to represent walking; and dynamic axial rotation (0.9°) for periodic movements during various daily activities of life (both at 0.5 Hz). Per activity, the distributions of average deviatoric shear strain εdev (defined as ) and fluid velocity vf were determined. A time average magnitude of each parameter per activity was used for a daily activity level, resulting in a MS ψ [23,36] (Eq 2)
(Eq 2)
where the constants, i.e. 5.5 for standing, 2 for walking, 4.5 for periodic movements and 4 for sitting resp., represent the time weighted average per activity in hours. Finally, to determine the daily activity level, the summation of εdev and vf was divided by 16 hours. The MS determines the preferred phenotype per integration point. The theory of Prendergast et al. [23,36] was adjusted to account for baseline residual strains in the unloaded disc (swelling-collagen tension equilibrium) [22]. This was based on the assumption that the original mechanoregulation theory did not take residual strains into account and that these strains are on the same order of magnitude as load-induced strains in IVD tissues where substantial osmotic swelling is balanced by collagen tensile strains. Initially, the threshold for strain and fluid velocity between the different preferred cell phenotypes was adjusted until the MS inside the geometrically defined nucleus region was mostly cartilage-favoured, i.e. the threshold between cartilage and fibrotic preferred phenotype was shifted in deviatoric strain from 11.25% to 14%; fluid velocity threshold was not adapted. Therefore, Eq 2 deviates from the original mechanoregulation theory [36], i.e.
instead of
.
Dependent on the preferred phenotype, cells went into proliferation, differentiation or into apoptosis and produced or degraded matrix, according to (Eq 3).
(Eq 3)
in which ci is the concentration of cell type i, i.e. chondrocytes or fibroblasts,
is the proliferation rate of cell type i, cspace for the maximum allowed cell concentration minus the current concentration,
for the differentiation rate of cell type i and
for the apoptosis rate of cell type i. The proliferation, differentiation and apoptosis rates for each cell type were constants and were either on or off, i.e. according to the preferred phenotype according to the MS. The different rates for the two cell types were previously calculated from experimental data collected in an extensive literature review [37] (Table 2).
Proliferation differentiation and apoptosis rate of fibroblast and chondrocytes chondrocytes were previously calculated in normalized fashion from experimental data collected in an extensive literature review [37].
Production and degradation of matrix was modeled as cell based. Both cell types could produce and degrade collagen and GAG although the rates at which this is possible were different between the two cells. The rate, at which matrix was produced, was dependent on the MS and the current cell concentration (Eq 4).
These equations contained for both cell types several components: 1) a synthesis part where i stands for GAG and collagen resp., in which the rate was dependent on cell type, with corresponding fixed rate (literature based), as well as the preferred target value, e.g. if the current amount of GAG is very small compared to the preferred amount, the synthesis rate is high; 2) two types of degradation: degrad1 is degradation of matrix by the non-preferred cell phenotype (literature based, Table 3), degrad2 is a degradation rate dependent on the current amount of matrix present compared to the preferred amount (fixed rate: 0.0011 (day -1), independent of cell phenotype. Since no literature value was available for the latter, it was assumed to be in the same range as degrad1 rate. If the current amount was much higher than the preferred concentration, the degrad2 rate was higher. For example, when going from fibrous tissue to cartilage, chondrocytes produce collagen and GAG and fibroblast degrade their matrix, i.e. collagen and GAG. However, since the concentration of collagen is too high for cartilage, chondrocytes will also break down collagen. In addition to this, there is a natural turn-over by chondrocytes and fibroblast which continuously produce and degrade their own matrix. For healthy IVD tissue, i.e. both healthy NP and AF tissue, the half-life time of collagen is approx. 95 years [38] (0.00003 day -1), the half-life time of GAG is 12 and 11.2 years for NP and AF tissue, respectively (0.00016 day -1) [39]. To account for this, Eqs 4 and 5 were extended with a continuous production and degradation of GAG and collagen by the preferred cell type.
The represents the preferred matrix content of the fibroblast;
represents the preferred matrix content of chondrocytes.
The process of tissue adaptation then was repeated several times until no more changes in disc composition were found, i.e. reaching the equilibrium state. It was assumed though that the orientation of the collagen fibres in the AF will not change due to adaptation. This equilibrium state was used as starting point to investigate the effect of surgical intervention. Due to surgical intervention the boundary conditions (BC) change, caused by a change in kinematic behaviour, resulting in a different MS per integration point (IP) and subsequently in a different ECM composition per IP (Fig 1). An explicit time integration scheme is implemented to simulate the changes in disc properties over time (time step of one day).
Patient-specific implementation.
An FE model of the patient disc was created based on the geometry obtained from the patient MR-scan. The composition of the disc was based on the Pfirrmann grade using the same MR-scan. Since the calculation of the equilibrium situation as described above is very time consuming and requires knowledge of patient-specific loading conditions, it was not possible to do this for each patient. To be able to predict patient-specific changes, an alternative approach was used. With this approach, the patient-specific disc model was used only to calculate the element composition: the initial solid volume fraction ns,0, the fibril density , and the fixed charge density FCD. These parameters were then mapped to the mechanoregulated disc with generalized geometry (Fig 2). This mapping is possible because all patient models are derived from a fixed mesh template that is morphed to the patient based on anatomical landmarks [44,45]. Since all elements thus have the same anatomical position, it is possible to transfer information from the patient-specific to the generalized model. The simulation of the mechanoregulation process then is performed using the generalised model, to predict changes in the disc composition, and the results are mapped back to the patient-specific model where they were used to calculate the composition-dependent material properties. The sequence of loading conditions representing a normal daily loading cycle was applied again, but scaled in magnitude according to the change in loading obtained from the lumbar spine model, before and after fusion. Directly after applying the adapted loading configuration due to fusion surgery to the model, the tissue changes were expected to be the largest. Therefore, the tissues were allowed to remodel for 3 months before the model was re-run to determine the new mechanical stimulus. After the first year, smaller changes were expected as the largest adaptation took place during the first year, the time step was set to 6 months.
To make predictions for a patient, first a patient-specific model is generated, including the patient specific composition, geometry and changes in loading conditions. These properties are transferred to the standard model using a mapping procedure. With this procedure, changes in any of the parameters are translated to equivalent changes in parameters used for the standard model. The results of the adaptation simulation are mapped back to the patient model.
Bone adaptation framework
To account for changes in the adjacent vertebra, i.e. changes in bone stiffness, the model was extended with a previously developed bone remodelling theory [28]. This theory is a homogenized version of an earlier microstructural bone remodelling theory [46,47] and predicts the bone density as a function of mechanical loading, quantified by the strain energy density (SED), and parameters related to cell-level activity and microstructure. Similar to the disc adaptation theory, the actual bone adaptation takes place at the level of the cells and is assumed to be regulated by the tissue-level SED at the surface of the bone tissue ST as sensed by the cells. Based on this mechanical signal, the osteocyte sensitivity μ and the osteoblast bone formation rate τ, the amount of newly formed bone is calculated. Also at this level, osteoclast cells were assumed to remove bone tissue at a constant speed determined by the osteoclast resorption chance focl and resorption volume Vres. Based on these parameters and the available bone surface per volume BS/TV, of which a fraction α was assumed to be available for bone remodelling, modelled as cylindrical structures, the bone resorption rate is calculated. By integrating these rates over the bone volume, the total change in bone density BV/TV (dimensionless, range 0–1.0) then is calculated as:
(Eq 5)
where ConstSED is defined as
(Eq 6)
As with the disc adaptation, an explicit time integration using fixed time step of one day is used to predict the changes in bone density over time.
Patient-specific implementation.
For the patient-specific implementation, the bone geometry and density distribution were based on the direct pre-operative patient CT scan. Whereas the bone remodelling simulation was performed for the full vertebrae, results will focus on changes in the bone density of the vertebral core. Similar as for the disc, applying the remodelling equations (Eqs 5 and 6) would typically require a time consuming tuning procedure. This could be solved by using a mapping procedure, similar to that used for the disc. A more elegant solution in this case, however, is to use the pre-operative state as the reference state and to take the changes before and after the operation as the driving force for remodelling (‘site specific‘ bone remodelling [25]) (Fig 3). The bone remodelling theory therefore was modified to enable the prediction of changes in bone density as they relate to changes in bone loading by defining the stimulus as the relative difference in the SED distribution after and before the fusion operation:
(Eq 7)
The pre-operative SED, or SEDref, distribution was calculated based on the direct pre-operative density distribution as obtained from the patient CT scan using loading conditions obtained from the full lumbar model (see below) representing the spine before the operation. The post-operative SED was calculated using the bone density distribution predicted by the remodelling theory using loading conditions obtained from the full lumbar model with fused disc. If both are similar, no remodelling will occur, i.e. Eq 7 and Eq 5 will be zero resp. Directly after applying the adapted loading configuration resulting from fusion surgery, the tissue changes were expected to be the largest. Therefore, the tissues were allowed to remodel for 3 months before the model was re-run to determine the new mechanical stimulus. After the first year, smaller changes were expected as the largest adaptation took place during the first year, the time step was set to 6 months. For the bone remodelling only the compressive load case was considered since this was found to be the most dominant loading mode and inclusion of other loading modes did not change the results much.
From CT scan of a patient, bone density and the geometry are obtained. From the full lumbar model patent-specific boundary conditions and changes in these boundary conditions due to fusion are obtained. The patient-specific FE model is used to determine the organ level load distribution. Using the bone remodelling theory, this is translated to a change in the density distribution.
In order to obtain realistic time-dependent bone remodelling patterns, some of the bone remodelling parameters were tuned. For this tuning, the bone density of one specific patient (patient #10, see below) was used as input. The remodelling process then was simulated using different values for the remodelling parameters until good agreement was achieved between the predicted 2-years results and the 2-years results obtained from the CT-scan. The values of the material parameters used is listed in Table 4.
Lumbar model
The lumbar model is an FE model representing the full lumbar part (L1 –S1) of the spine [29] (Fig 4). This model, which includes vertebrae, discs and ligament, is generated from patient CT and MR scans and has a patient-specific geometry of both the vertebrae and discs [44,45]. Each vertebra was meshed with around 40,000 and each disc with around 20,000 8-node hexahedron elements, such that the total number of elements was around 300,000. Typical element size was between 1.5 and 2.5 mm.
An example of a lumbar FE model. The model contains vertebrae, discs and ligaments form L1 till S1.
The cortical bone and bony endplates were represented by a structural mesh layer at the periphery of the vertebral body models. The thickness of this layer varied depending on the location, i.e. middle, cranial/caudal, anterior or posterior cortical bone, and middle or peripheral bony endplate anterior posterior, and local thickness values were defined from direct measurements on histological cuts [48]. The material properties were chosen linear elastic, with a Young's modulus of 12 GPa in the axial direction and 8 GP in the transversal direction for the cortical bone and a Young's modulus of 250 MPa in the axial direction and 140 MPa in the transversal direction for the cancellous bone [16]. The AF and the NP were both modelled as poro-hyperelastic materials [29]; different to what is used for disc adaptation. The total stress tensor σ caused by external loads was the superimposition of the porous solid stress and the fluid pore pressure, that were respectively derived from a strain energy density function W, and Darcy's law [29].
(Eq 8)
(Eq 9)
(Eq 10)
where G and K are the shear and bulk modulus respectively of the porous solid (drained), J the determinant of the deformation tensor F, I1 the first strain invariant, I is the second order unit tensor, uf the pore fluid velocity and ϕ and k are the porosity and hydraulic permeability tensor resp. The term WANI is an anisotropic strain energy density term, different from zero only for the AF. For more details, see Malandrino et al. [29].
The model included the seven major ligaments (anterior longitudinal, posterior longitudinal, capsular, supraspinous, ligamentum flavum, interspinous and intertransverse ligaments) as described before [22]. In summary, the ligaments, except for the posterior longitudinal ligament (PLL), were modelled with truss elements with an exponential stress-strain relationship for the toe region, followed by a linear behavior [16] as described by
where S is the 2nd Piola-Kirchhoff stress, E the Green strain and A, B, C and D four constants defined in Noailly et al. [49]. The PLL is not only attached to the vertebrae but also to the annulus. Thus, it was modified to match the AF posterior shape, and modelled with surface fibre-reinforced elements, tied to the edge of the vertebrae and the posterior annulus. The material properties of the rebars, i.e. the reinforcing fibres, were hypoelastic similar to the other ligaments. The articular facet surfaces were modelled with surface-to-surface frictionless contact. A similar model of a similar motion segment, but with a generic geometry, was validated earlier relative to experimental measurements and provided a good representation of the time-dependent decrease of the reaction forces, the range of motion for applied moments and facet forces [22].
Load boundary conditions, i.e. compressive forces and flexion moments, based on the patients’ weight and height, respectively, were simulated. Resultant forces (follower loads) and forward bending moments were applied at the centre of the cranial endplates, according to the musculoskeletal model calculations and anthropometry-dependent interpolations reported by Han et al [50]. In flexion, the effects of the imposed moment were always analyzed at 12o of calculated rotation. The lumbar spine model was subsequently used to investigate the resulting load conditions per segment with and without simulated fusion. These local loads were used to calculate the percentage change in loading for the segment adjacent to fusion, which determined the proportion by which the daily loading cycle should be adjusted for the disc adaptation model, so as to reflect the mechanical effect of the fusion simulated patient-specifically.
Integrated bone–disc adaptation
To study the effect of fusion surgery on both bone and disc simultaneously, the two algorithms were integrated into the same framework. Bone remodelling is assumed to be driven by the magnitude of the applied load [51], i.e. the magnitude of a peak load, while disc adaptation is assumed to be dependent on the average load distribution during normal daily activity (see disc adaptation framework description).
Since bone remodelling occurs much faster than disc degeneration (order of months versus order of years), the simulation of bone remodelling requires much smaller time steps than that of disc adaptation. Using a fully coupled approach would require the evaluation of both the stresses/strains in the disc and in the bone at each time step. Since evaluating the average stresses/strains in the disc is very computationally expensive, and since the changes in disc composition would be minor for such small time steps, a semi-coupled approach was used instead. With this approach, the two processes where separated and performed in series to each other once for each larger time step, i.e. first bone remodelling followed by disc adaptation (Fig 5). Within each larger time step, bone remodelling simulations were allowed to come to their temporary steady state; as disc adaptation is a slower process, the adaptation was simulated in a fixed time frame not necessarily resulting in steady state. Before tissue adaptation was simulated (both disc and bone), the disc was allowed to come to mechanical steady-state by applying three preconditioning loading cycles consisting of an axial compression load [31,32,33,34,35] followed by a creep period, similar as described in the disc adaptation framework. This precondition mimics the action of external loads as well as muscle activation.
In the coupled model, the bone remodelling and the disc adaptation algorithms are combined to represent the bone remodelling and disc adaptation due to intervention. Starting with the bone remodelling, a new density distribution is calculated that is mapped to the disc adaptation model. Changes in disc composition are calculated and are subsequently mapped back to the patient-specific model. Each module is explained in detail above.
Patient study
To demonstrate the feasibility of this described generic framework to investigate changes in adjacent tissue after a fusion, it was applied to the data of 10 patients (pt) suffering from low back pain. The original database was built up in the EU-funded MySpine project and contained the prospectively collected clinical and imaging data of 192 patients treated at the National Center for Spinal Disorders according to the national guidelines and institutional protocols. This retrospective study was based on a subset of the database and approved by the Scientific and Research Ethics Committee of the Hungarian Medical Research Council (751/PI/2010). Patients with known metabolic bone disease were excluded from the study cohort. This cohort consisted of young adults (mean age: 38 years) where osteoporosis is rare. The cohort did not contain patients with clinically poor bone quality. Based on CT and MRI data, pt-specific FE meshes of the lumbar spine were reconstructed [29,44,45]. Out of these 192 cases, 27 patient cases underwent a monosegmental lumbar fusion. For 10 cases (7/10 male; mean age at surgery 38 years; 5/10 fusion at L4-L5, 4/10 fusion at L5-S1, 1/10 L3-L4), the lumbar spine model [29] could successfully solve all applied loading configurations (axial compression and flexion for the pre-surgical situation and after fusion treatment resp.). For these cases, the effect of fusion surgery on the adjacent disc and vertebra was simulated and compared to the clinical follow-up data (Fig 6).
From the lumbar model the two vertebrae neighbouring the degenerated disc and the neighbouring disc are extracted as the ‘patient-specific’ model. Also the boundary conditions representing the pre-operative loading are extracted from the lumbar model. Following, the fusion surgery was simulated and new boundary conditions were calculated. For these new boundary conditions, changes in the disc and the vertebrae neighbouring the operated region were calculated.
Outcome measures
Changes in bone density of the vertebra and changes in ECM at the adjacent level were compared to clinical follow-up (FU) data 2 years post-surgery. To compare the predicted bone densities, the simulated bone density value per location was compared to the clinically observed bone density at the same location 2 years post-surgery. This was done for the vertebra adjacent to the fusion level, e.g. when level L4-L5 was fused, the predicted and clinical densities of vertebra L3 were compared to each other. To do so, the same mapping and morphing procedure as to determine the initial geometry and bone density [44,45] was used for the 2 year FU data. To reduce the effects of noise in the CT scan and to obtain more meaningful bone density values, the element bone density was evaluated for a spherical region with a radius of 2 mm with its centre at the element centre [52]. The calculated density then was assigned to all integration points of the element. As elements were typically smaller than 2 mm, this resulted in overlapping evaluation regions when calculating densities of neighbouring elements, which naturally smoothens the density distribution. This was done for both simulated and clinical FU data. Since the reproducibility of the results depends on the reproducibility of the meshing algorithm, a small sub-study was performed to investigate errors that can be expected due to meshing errors. In this sub-study, meshes for four patients were generated three times, and the density distribution was calculated. Based on these results it was concluded that the mean detectable change (MDC) per element was 1.9%, indicating that changes less than this value cannot be detected. At the periosteal boundary of the vertebrae, however, a much lower reproducibility was found. This was due to the fact that the meshes were generated by segmenting the vertebrae at a rather low density, thus to be sure all bone was meshed. As a consequence of this, however, elements at the periosteal side sometimes were just outside of the bone region, showing a very low density. This did not affect the results of the stiffness or the remodelling though. The similarity between the simulation and FU data was expressed by a correlation coefficient, i.e. the bone density in each element at the adjacent level of the fusion was compared to the density in the same element in the simulation. In order to have a clinically relevant outcome for the disc changes, the relative changes in ECM content were translated into a change in pt-specific degeneration grade for the overall IVD based on tissue water content. The relative change in ECM content between different degeneration grades has been well studied in literature [53,54] and the best predictor of a change in degeneration grade is water content within the tissue [53]. As in the IVD model the water content is based on the equilibrium between collagen tension and GAG content-based osmotic swelling, changes in collagen and GAG content in IVD tissue lead to a change in water content. Due to the used modelling approach this predicted change in water content is a relative change. Prior to treatment, for each patient, the grade of degeneration was known, and thus the water content inside the IVD tissue [53]. After adaptive simulation, the calculated relative change in water content was translated to an absolute water content, based on the prior known pt-specific water content and the simulated change (grade dependent, based on literature [53], leading to the pt-specific change in water content and corresponding pt-specific change in grade. This was compared to clinical scored degeneration grade per patient 2 years post-surgery. Three randomly selected cases (#4, #6 and #9) were simulated up till 10 years post-surgery to study the effect of fusion on both adjacent vertebra and IVD at a longer time span. Although these results cannot be compared to clinical data, it shows the potential of the developed framework.
Results
Disc steady state
To verify the adaptation method and to determine a steady-state IVD with corresponding matrix content and cell density spatial distribution, it was assumed that the adaptation changes for physiological loads should result in a steady-state IVD with tissue morphologies of a healthy normal IVD. After adaptation to steady state, in the centre of the IVD, a large region was present where cartilage tissue had formed. Circumferentially in the outer regions of the IVD, steady-state fibrous tissue had formed with corresponding matrix components. In between these two regions, particularly at posterior location of the IVD, a mixed tissue is present (Fig 7).
Cross-sections of a steady state disc with (a) collagen (b) GAG and (c) water (relative amount of content). In the centre of the disc, a region with low collagen (blue) and high GAG (red) and water content is found, similar to cartilage tissue, i.e. NP. Towards the outer shape of the disc, regions with high collagen and low GAG and water content are observed, i.e. AF tissue. In between these regions, a mixture of both tissues is observed, indicating the transition zone between both tissues. The black lines indicate the borders between the original geometrically defined NP and AF region.
Pt-specific simulations
Bone adaptation.
In general, modest changes in the adjacent tissues were found, both in bone and IVD tissue. For bone, the predicted bone densities were in fair to good agreement with the clinical FU data (Table 5, Fig 8). The low density cortical appearance of the vertebrae in Fig 8 relates to the fact that elements at the most periosteal locations can be outside of the vertebral body. As mentioned earlier, however, this will not affect the results. The fact that these elements in some cases were just outside and some other case just inside the bone, however, did considerably reduce the correlations as calculated in Table 5. Note, in this section, the visual results of bone remodelling are visualized on the pt-specific adjacent vertebra only, i.e. the caudal vertebra of the segment adjacent to the fusion without disc.
Examples of graphical representation of bone density (BV/TV) as predicted 2 years post-surgery (b, e and h) and measured at the same time point (c, f and i), compared to pre surgery (a, d and g). a, b and c correspond to patient #1, d, e and f to patient #7, g, h and i to patient #10. A black colour indicates a density of 0. Note that the dark-blue appearance at the cortical shell is because the mesh boundaries were taken slightly larger than the vertebral size. As a result, the outer layer of elements can have a very low density.
Correlation between predicted bone density (BV/TV) vs clinical FU bone density. Case #10 was used to fit the remodelling parameters, resulting in a very high correlation.
Disc adaptation.
As for the adjacent disc, in all patients, only modest changes in biochemical content were predicted by the simulation. When looking at the relative amount of collagen, the nucleus region had become more fibrous, indicating that degeneration was progressing (see Fig 9 for typical results, note that to clearly show the effect of change in loading, the results are visualized on the generic disc geometry). Similar changes but in the reverse direction was observed for the GAG content as well as the water content. Although the changes were modest, the results indicated that the biochemical content was progressing towards the next stage of degeneration. However, based on the relative change in water content inside the IVD, no grade change would be anticipated in the upper level adjacent to the treated disc within the first 2 years after surgery. Matching this result, no change in degeneration grade was observed clinically in all cases.
a) Start collagen content and predicted collagen content 2 years after fusion (b). c) and d) Change in GAG content for 2 year simulation due to fusion. e) and f) Changes in water content in the disc for 2 year simulation. For visualization purposes, only the disc is shown. Predicted changes are for case #9 and serve as example of the outcome for the other 9 cases.
Three cases (#4, #6, #9) were simulated up till 10 years. For the first two years, the predicted tissue changes correspond well to the clinical data and only modest changes were predicted (Table 5). For case #9, after 5 years a significant change in ECM content was simulated as well as a significant change in bone density (Fig 10, Fig 11 and Fig 12). For the other two cases, only modest changes in both bone and ECM content were simulated. After 10 years, in all cases an increase in degeneration of one grade was predicted.
Tissue changes in both bone and disc changes simulated up till 10 years post-operative in a combined view. For disc tissue, the ECM components are visualized: collagen, GAG and water content (relative amount). Water content is visualized in black and white, mimicking the water content as could be observed on MRI data. Next to these ECM changes, also bone changes are visualized for the 4 different time points. After 2 years, the clinical FU-data is visualized to show a graphical comparison. Bone density values vary between 0 and 1, relative collagen content between 0.65 and 0.15, GAG between 0.85 and 0.35, water content between 0.81 and 0.75 respectively.
Tissue changes in both bone and disc changes simulated up till 10 years post-operative in a combined view. For disc tissue, the ECM components are visualized: collagen, GAG and water content (relative amount). Water content is visualized in black and white, mimicking the water content as could be observed on MRI data. Next to these ECM changes, also bone changes are visualized for the 4 different time points. After 2 years, the clinical FU-data is visualized to show a graphical comparison. Bone density values vary between 0 and 1, relative collagen content between 0.65 and 0.15, GAG between 0.85 and 0.35, water content between 0.81 and 0.75 respectively.
Tissue changes in both bone and disc changes simulated up till 10 years post-operative in a combined view. For disc tissue, the ECM components are visualized: collagen, GAG and water content (relative amount). Water content is visualized in black and white, mimicking the water content as could be observed on MRI data. Next to these ECM changes, also bone changes are visualized for the 4 different time points. After 2 years, the clinical FU-data is visualized to show a graphical comparison. Bone density values vary between 0 and 1, relative collagen content between 0.65 and 0.15, GAG between 0.85 and 0.35, water content between 0.81 and 0.75 respectively.
Discussion
The constituent based FE model of an IVD was successfully extended to a mechanoregulated IVD model. Based on a mechanoregulated tissue differentiation theory [23], a stable healthy IVD with a realistic ECM spatial distribution was obtained. Subsequently, the IVD model was incorporated into a spinal segment model and was extended with a bone remodelling theory [28]. Finally, this method was successfully extended to a patient-specific algorithm and demonstrated for 10 cases that underwent fusion surgery to predict vertebral bone density and disc tissue adaptive changes in the adjacent segment.
In both the clinical study and the computer simulations, a modest increase in vertebral density was found as the result of a neighbouring fusion operation. As observed in multiple studies, this effect may lead to accelerated disc degeneration of the adjacent disc levels since an increase in stiffness may lead to changes in loading configuration at the disc as suggested by Pye et al. [13]. In the current study, however, the observed and simulated changes in density were small, as only a 2 year follow-up period was available. Interestingly, the simulation study predicted ECM changes towards the next degeneration grade, suggesting that a longer follow-up period would reveal more severe degeneration. This was clear also from the three cases in which simulation results were continued until 10 years follow-up, where a clear trend in disc changes were observed, especially from 5 years post-surgery, even when changes within the first 2 years were modest. These findings are in line with those found by Pye et al. [13], Aota et al. [8] and Etebar and Cahill [55] who observed a change in grade after 25 months or more.
A point that requires further discussion is the fact that the total adaptation model as presented here involves a large number of parameters, of which only a limited number (i.e. disc/vertebral geometry, disc composition, bone density distribution and changes in load between the pre- and post-operative situation) could be determined in a patient-specific manner. The remaining ones were based on literature values or determined using a fitting procedure. The solution of this fitting procedure may not be unique. For example, it was shown in earlier studies that the results for the disc are not sensitive for all parameters [22], making them difficult to fit. Similarly, for the bone remodelling simulations, different combinations of the osteocyte sensitivity μ and the osteoblast bone formation rate τ can provide the same results, and due to the specific formulation used the results are not dependent on the osteoclast resorption chance focl and resorption volume Vres (when assuming these do not change pre- and post-operative). Although it would be possible to remove or lump such parameters, it was chosen not to do so here since they all have a clear physical interpretation such that reasonable estimates are possible from literature. Keeping them included may make the model more versatile in case of future improvements in patient measurements, and enable the model to account for changes that are neglected in this study (e.g. effects of bisphosphonate treatment that affect osteoclast activity). Whereas the sensitivity of the results of the disc degeneration and the bone remodelling algorithm has been investigated in earlier studies [22,28], it was not possible to address this for the combined patient-specific model used here. As such, the present study merely serves as a demonstration of the feasibility of such a combined model than as a rigorous validation. Obvious, a parameter sensitivity study would be required to better investigate the effect that individual parameters have on the outcome and thus the need for patient-specificity. The bone remodelling algorithm taken from Colloca et al., [28] was modified to account for the difference in SED values in the pre- and post-operative situation. The actual remodelling algorithm used here is similar to those used to predict bone loss around implants (e.g. Huiskes et al., [25]), and also referred to as a ‘site-specific’ remodelling. The major advantage of this approach is that realistic results can be obtained even when evaluating only a limited number of daily loading conditions and that it results in more robust predictions.
Although predicted changes in bone density show similar trends as those observed in the patients, there also are clear differences. Several reasons could explain such differences. First, the change in loading due to fusion was calculated from the full lumbar model [29] only once by comparing the pre- operative and post-operative loads. Ideally, the BCs are re-determined after each adaptation cycle since changes in disc adaptation might affect spine behaviour. However, such effects are expected to be very small. Second, the value of some remodelling parameters was determined by fitting the simulation results and the clinical results of only one patient. A more rigorous fitting involving more patients and a more advanced optimization scheme might improve the overall performance of the model. Since the present study involved only 10 patients, this was not possible. Third, a standardized loading protocol was used for all patients and it was assumed that the activity levels before and after the operation would be the same. Although the model could account for different loading patterns, e.g. related to occupational behaviour, or changes in activity level after the operation, no information was available for the patients included here. Fourth, there might be patient-specific differences in bone and disc metabolism and adaptation rates. For example, in patients on bisphosphonate treatment the changes in bone density might be less than expected. Since the model includes the cell level, it could also account for such changes in bone and disc responsiveness, provided that the effects of drugs on cell activity are known. For the present study, no information about factors that may affect the bone or disc adaptation was available.
There were also a number of technical limitations to the model used here. Most importantly, it was not possible to measure patient-specific disc composition based on the MR images. Consequently, only the Pfirrmann grade could be used to estimate the disc composition. The present study, however, was not powered to detect significant changes in Pfirrmann grade within a 2 years follow-up period. As a result, no changes in Pfirrmann grades could detected and thus it was not possible to validate the disc degeneration model relative to the patient study. As mentioned above, this limits the present study to a feasibility study. Possibly, more advanced MR techniques, such as T1rho and T2 mapping may better reveal the patient-specific disc composition.
Second, because of the lack of accurate patient-specific disc composition and because the calculation of the equilibrium situation for the disc is very time consuming, it was not possible to get an equilibrium state for the disc for each individual patient. Instead the equilibrium situation was determined only for a generalized disc geometry and a mapping procedure was used to map predicted compositional changes back to the patient model. Since the generalized model obviously lacks the specific patient geometry of disc and vertebrae, this approach cannot account for geometry-related (changes in) force distribution. It does, however, account for composition-dependent (changes in) load transfer.
Finally, in order to predict a realistic material distribution, some parameters in the mechanoregulation theory had to be modified. This was based on the assumption that the original theory did not take residual strains into account and that these strains are on the same order of magnitude as load-induced strains in IVD tissues where substantial osmotic swelling is balanced by collagen tensile strains. The accuracy of the parameters as found here is not known since no other data is available to which these can be compared. As such, they should be treated as 'effective' parameters only.
Conclusion
In conclusion, a framework was developed and tested that predicts patient specific tissue changes adjacent to the fusion region. To the best of our knowledge, this is the first framework able to simulate pt-specific tissue changes after lumbar spinal fusion in both vertebra and disc. Both bone and disc tissue changes correlated with the clinical observed findings at 2 years post-surgery.
References
- 1. Andersson GBJ (1999) Epidemiological features of chronic low-bak pain. The Lancet 354:581–585
- 2. Hoy D, Brooks P, Blyth F, Buchbinder R (2010) The epidemiology of low bak pain. Best Practice & Research Clinical Rheumatology 24:769–781
- 3. Cheung KM, Karppinen J, Chan D, Ho DW, Song YQ, Sham P, et al. (2009) Prevalence and pattern of lumbar magnetic resonance imaging changes in a population study of one thousand forty-three individuals. Spine 34:934–940 pmid:19532001
- 4. Luoma K, Riihimaki H, Luukkonen R, Raininko R, Viikari-Juntura E, Lamminen A (2000) Low back pain in relation to lumbar disc degeneration. Spine 25:487–492 pmid:10707396
- 5. Airaksinen O, Brox JI, Cedraschi C, Hildebrandt J, Klaber-Moffett J, Kovacs F, et al. (2006) European guidelines for the management of chronic non-specific low back pain. European Spine Journal 15:192–300
- 6. Carreon LY, Glassman SD, Djurasovic M, Dimar JR, Johnson JR, Puno RM, et al. (2009) Are preoperative health-related quality of life scores predictive of clinical outcomes after lumbar fusion? Spine 34:725–730 pmid:19333106
- 7. Harrop JS, Youssef JA, Maltenfort M, Vorwald P, Vaccaro AR, Hilibrand AS (2008) Lumbar adjacent segment degeneration and disease after arthrodesis and total disc arthroplasty. Spine 33:1701–1707 pmid:18594464
- 8. Aota Y, Kumano K, Hirabayashi S (1995) Postfusion instability at the adjacent segments after rigid pedicle screw fixation for degenerative lumbar spinal discorders. Journal of Spinial Disorders 8:464–473
- 9. Kumar MN, Baklanov A, Chopin D (2001) Correlation between sagittal plane changes and adjacent segment degeneration following lumbar spine fusion. European Spine Journal 10:314–319 pmid:11563617
- 10. Lee CK (1988) Accelerated degeneration of the segment adjacent to a lumbar fusion. Spine 13:375–377 pmid:3388124
- 11. Radcliff KE, Kepler ChK, Jakoi A, Sidhu GS, Rihn J, Vaccaro AR, et al. (2013) Adjacent segment disease in the lumbar spine following different treatment interventions. The Spine Journal 13:1339–1349 pmid:23773433
- 12. Adams MA, Pollintine P, Tobias JH, Wakley GK, Dolan P (2006) Intervertebral disc degeneration can predispose to anterior vertebral fractures in the thoracolumbar spine. Journal of Bone and Mineral Research 21:1409–1416 pmid:16939399
- 13. Pye SR, Reid DM, Adams JE, Silman AJ, O'Neill TW (2006) Radiographic features of lumbar disc degeneration and bone mineral density in men and women. Annals of the Rheumatic Diseases 65:234–238 pmid:16014671
- 14. Wang Y, Boyd SK, Battie MC, Yasui Y, Videman T (2011) Is greater lumbar vertebral BMD associated with more disk degeneration? A study using uCT and Discography. Journal of Bone and Mineral Research 26:-2785 pmid:21786320
- 15. Malandrino A, Planell JA, Lacroix D (2009) Statistical factorial analysis on the poroelastic material properties sensitivity of the lumbar intervertebral disc under compression, flexion and axial rotation. Journal of Biomechanics 42:2780–2788 pmid:19796766
- 16. Noailly J, Wilke H-J, Planell JA, Lacroix D (2007) How does the geometry affect the internal biomechanics of a lumbar spine bi-segment finite element model? Consequences on the validation process. Journal of Biomechanics 40:2414–2415 pmid:17257603
- 17. Ruberté LM, Natarajan RN, Andersson GBJ (2009) Influence of single-level lumbar degenerative disc disease on the behavior of the adjacent segments: a finite element model study. Journal of Biomechanics 42:341–348 pmid:19136113
- 18. Schmidt H, Heuer F, Simon U, Kettler A, Rohlmann A, Claes L, et al. (2006) Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus. Clinical Biomechanics 21:337–344 pmid:16439042
- 19. Schmidt H, Shirazi-Adl A, Galbusera F, Wilke H-J (2010) Response analysis of the lumbar spine during regular daily activities—a finite element analysis. Journal of Biomechanics 43:1849–1856 pmid:20394933
- 20. Homminga J, Aquarius R, Bulsink VE, Jansen ChTJ, Verdonschot N (2012) Can vertebral density changes be explained by intervertebral disc degeneration. Medical Engineering & Physics 34:453–458
- 21. van Rijsbergen MM, Barthelemy VMP, Vrancken ACT, Crijns SPM, Wilke H-J, van Rietbergen B, et al. (2016) Moderately degenerated lumbar motion segments: are they truly unstable? Biomechanics and Modeling in Mechanobiology
- 22. Barthelemy VMP, van Rijsbergen MM, Wilson W, Huyghe JM, van Rietbergen B, Ito K (2016) A computational spinal motion segment model incorporating a matrix composition based model of the intervertebral disc. Journal of the Mechanical Behavior of Biomedical Materials 54:194–204 pmid:26469631
- 23. Prendergast PJ, Huiskens R, Soballe K (1997) Biophysical stimuli on cells during tissue differentiation at implant interfaces. Journal of Biomechanics 30:539–548 pmid:9165386
- 24. Isaksson H, Comas O, van Donkelaar CC, Mediavilla J, Wilson W, Huiskens R, et al. (2007) Bone regeneration during distraction osteogenesis: Mechan-regulation by shear strain and fluid velocity. Journal of Biomechanics 40:2002–2011 pmid:17112532
- 25. Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Sloof TJ (1987) Adaptive bone-remodeling theory applied to prosthetic-design analysis. Journal of Biomechanics 20:1135–1150 pmid:3429459
- 26. Huiskes R, Weinans H, Dalstra M (1989) Adaptive bone remodeling and biomechanical design considerations for noncemented total hip arthroplasty. Orthopedics 12:1255–1267 pmid:2798252
- 27. Van Lenthe GH, de Waal Malefijt MC, Huiskes R. (1997) Stress shielding after total knee replacement may cause bone resorption in the distal femur. J Bone Joint Surg Br. Jan;79(1):117–22 pmid:9020459
- 28. Colloca M, Ito K, van Rietbergen B (2014) An analytical approach to investigate the evolution of bone volume fraction in bone remodeling simulation at the tissue and cell level. Journal of Biomechanical Engineering 136:031004-1-031004-8
- 29. Malandrino A, Pozo JM, Castro-Mateos I, Frangi AF, van Rijsbergen MM, Ito K, et al. (2015) On the relative relevance of subject-specific geometries and degeneration-specific mechanical properties for the study of cell death in human intervertebral disk models. Frontiers in Bioengineering and biotechnology 3:1–15
- 30. Narmoneva DA, Wang JY, Setton LA (1999) Nonuniform swelling-induced residual strains in articular cartilage. Journal of Biomechanics 32:401–408 pmid:10213030
- 31. Adams MA, Dolan P, Hutton WC, Porter RW (1990) Diurnal changes in spinal mechanics and their clinical significance. Journal of Bone and Joint Surgery 72:266–270 pmid:2138156
- 32. McMillan DW, Garbutt G, Adams MA (1996) Effect of sustained loading on the water content of intervertebral discs: implications for disc metabolism. Annals of the Rheumatic Diseases 55:880–887 pmid:9014581
- 33. Paajanen H, Lehto I, Alanen A, Erkintalo M, Komu M (1994) Diurnal fluid changes of lumbar discs measured indirectly by magnetic resonance imaging. Journal of Orthopaedic Research 12:509–514 pmid:8064481
- 34. Rohlmann A, Graichen F, Weber U, Bergmann G (2001) Monitoring In Vivo implant loads with a telemeterized internal spinal fixation device. Spine 25:2981–2986
- 35. Wilke H-J, Neef P, Caimi M, Hoogland Th, Claes LE (1999) New In Vivo Measurements of pressures in the intervertebral disc in daily life. Spine 24:755–762 pmid:10222525
- 36. Huiskes R, van Driel WD (1997) A biomechanical regulatory model for periprosthetic fibrous-tissue differentiation. Journal of Materials Science: Materials in Medicine 8:785–788 pmid:15348791
- 37. Isaksson H, van Donkelaar CC, Huiskes R, Ito K (2008) A mechano-regulatory bone-healing model incorporating cell-phenotype specific activity. Journal of Theoretical Biology 252:230–246 pmid:18353374
- 38. Sivan SS, Wachtel E, Tsitron E, Sakkee N, van der Ham F, DeGroot J, et al. (2008) Collagen turnover in normal and degenerate human intervertebral discs as determined by the racemization of aspartic acid. The Journal of Biological Chemistry 283:8796–8801 pmid:18250164
- 39. Sivan SS, Tsitron E, Wachtel E, Roughley PJ, Sakkee N, van der Ham F, et al. (2006) Aggrecan turnover in human intervertebral disc as determined by the racemization of aspartic acid. The Journal of Biological Chemistry 281:13009–13014 pmid:16537531
- 40. Dimicco MA, Sah RL (2003) Depedence of cartilage matrix composition on biosynthesis, diffusion and reaction. Transport in porous media 50:57–73
- 41. Pisu M, Lai N, Concas A, Cao G (2006) A novel simulation model for engineered cartilage growth in static systems. Tissue Engineering 12:2311–2320 pmid:16968171
- 42. Sengers BG, van Donkelaar CC, Oomens CWJ, Baaijens FP (2004) The local matrix distribution and the functional development of tissue engineerd cartilage, a finite element study. Annals of Biomedical Engineering 32:1718–1727 pmid:15675683
- 43. Wilson ChG Bonassar LJ, Kohles SS (2002) Modeling the dynamic composition of engineered cartilage. Archives of Biochemistry and Biophysics 408:246–254 pmid:12464278
- 44. Castro-Mateos I, Pozo JM, Eltes PE, Del rio L, Lazary A, Frangi AF (2014) 3D segmentation of annulus fibrosus and nucleus pulposus from T2-weighted magnetic resonance images. Physics in Medicine & Biology 59:7847–7864
- 45. Castro-Mateos I, Pozo JM, Lazary A, Frangi AF (2015) 3D vertebra segmentation by Feature Selection Active Shape Model. Lecture notes in Computational vision and biomechanics 20:241–245
- 46. Ruimerman R, Hilbers P, van Rietbergen B, Huiskes R (2005) A theoretical framework for strain-related trabecular bone maintenance and adaptation. Journal of Biomechanics 38:931–941 pmid:15713314
- 47. van Oers RFM, Ruimerman R, Tanck E, Hilbers PAJ, Huiskes R (2008) A unified theory for osteonal and hemi-osteonal remodeling. Bone 42:250–259 pmid:18063436
- 48.
Noailly J, (2009) Model developments for in silico studies of the lumbar spine biomechanics, PhD Thesis. ed, Model developments for in silico studies of the lumbar spine biomechanics—PhD Thesis. Technical University of Catalonia, Barcelona, Spain
- 49. Noailly J, Ambrosio L, Elizabeth Tanner K, Planell JA, Lacroix D (2012) In silico evaluation of a new composite disc substitute with a L3-L5 lumbar spine finite element model. Eur Spine J. Jun;21 Suppl 5:S675–87
- 50. Han KS, Rohlmann A, Zander T, Taylor WR. (2013) Lumbar spinal loads vary with body height and weight. Med Eng Phys. 2013 Jul;35(7):969–77 pmid:23040051
- 51. Rubin CT, Lanyon LE (1985) Regulation of bone mass by mechanical strain magnitude. Calcied Tissue International 37:411–417
- 52. Harrigan TP, Jasty M, Mann RW, Harris WH (1988) Limitations of the contiuum assumption in cancellous bone. Journal of Biomechanics 21:269–275 pmid:3384825
- 53. Antoniou J, Steffen T, Nelson F, Winterbottom N, Hollander AP, Poole RA, Aebi M, Alini M (1996) The human lumbar intervertebral disc: evidence for changes in the biosynthesis and denaturation of the extracellular matrix with growth, maturation, ageing, and degeneration. The Journal of Clinical Investigation 98:996–1003 pmid:8770872
- 54. Johannessen W, Elliott DM (2005) Effects of degeneration on the biophasic material properties of human nucleus pulposus in confined compression. Spine 30:724–729
- 55. Etebar S, Cahill DW (1999) Risk factors for adjacent-segment failure following lumbar fixation with rigid instrumentation for degenerative instability. Journal of Neurosurgery: Spine 90:163–169