In scoliosis, kypholordos and wedge properties of the vertebrae should be involved in determining how stress is distributed in the vertebral column. The impact is logically expected to be maximal at the apex.
To introduce an algorithm for constructing artificial geometric models of the vertebral column from DICOM stacks, with the ultimate aim to obtain a formalized way to create simplistic models, which enhance and focus on wedge properties and relative tilting.
Our procedure requires parameter extraction from DICOM image-stacks (with PACS,IDS-7), mechanical FEM-modelling (with Matlab and Comsol). As a test implementation, models were constructed for five patients with thoracal idiopathic scoliosis with varying apex rotation. For a selection of load states, we calculated a response variable which is based upon distortion energy.
For the test implementation, pairwise t-tests show that our response variable is non-trivial and that it is chiefly sensitive to the transversal stresses (transversal stresses where of main interest to us, as opposed to the case of additional shear stresses, due to the lack of explicit surrounding tissue and ligaments in our model). Also, a pairwise t-test did not show a difference (n = 25, p-value≈0.084) between the cases of isotropic and orthotropic material modeling.
A step-by-step description is given for a procedure of constructing artificial geometric models from chest CT DICOM-stacks, such that the models are appropriate for semi-global stress-analysis, where the focus is on the wedge properties and relative tilting. The method is inappropriate for analyses where the local roughness and irregularities of surfaces are wanted features. A test application hints that one particular load state possibly has a high correlation to a certain response variable (based upon distortion energy distribution on a surface of the apex), however, the number of patients is too small to draw any statistical conclusions.
This paper concerns idiopathic scoliosis in adolescents. About 80 percent of structural scoliosis is classified as idiopathic scoliosis (American Association of Neurological Surgeons [1S. Aaro, and M. Dahlborn, "Estimation of vertebral rotation and the spinal and rib cage deformity in scoliosis by computer tomography", Spine, vol. 6, no. 5, pp. 460-467.[http://dx.doi.org/10.1097/00007632-198109000-00007] [PMID: 7302680] ]). If there is an underlying known diagnosis, believed to be the chief causing factor, then the scoliosis is usually not classified as idiopathic, an example of this is neuromuscular scoliosis. The pathological mechanisms of scoliosis are multi-factorial and obviously in the idiopathic case not yet fully understood. Main factors that are recognized as essential are rotational lordosis with a growth-discrepancy between the anterior and posterior, asymmetrical growth and muscle activities, see Hefti [2American Association of Neurological Surgeons, http:// www.aans.org/ Patient%20Information/ Conditions%20and%20Treatments/ Scoliosis.aspx]. Brink et al. [3M. Anburajan, and V. Divya, "Finite Element Analysis of Human Lumbar Spine", 2011 International Conference on Electronics Computer Technology, .] proposed that the anterior overgrowth is the result of the scoliotic mechanism. For deeper current knowledge about idiopathic scoliosis see e.g. the textbook by Newton et al. [4J.E. Akin, Finite Elements for Analysis and Design., Academic Press, .]. Using FEM-analysis based on models obtained from medical imaging, allows for noninvasive investigation of how stress is distributed at and near the apex, in scoliosis patients. Mechanical loading influences the rate of long bones and vertebrae and clearly altered loading is implicated in scoliosis (and other skeletal deformities), see e.g. Stokes [5H. Aleti, and I.A. Motaleb, "Study of Stress and Strain in Human Spine", IEEE International Conference, pp. 184-9.[http://dx.doi.org/10.1109/EIT.2014.6871759] ]. This together with the aforementioned observations of growth-discrepancy between the anterior and posterior in scoliosis makes it interesting to gather knowledge about how stress from applied axial and shear forces respectively, is distributed in the vertebral column at, and near the apex, for varied severity of scoliosis.
To find a semi-global approach to simplistic, but rigorously defined, modelling that keep the main wedge and relative tilt features of the vertebrae. The requirement of rigour and simplicity together implies, in practice, that small local surface variations will be discarded and instead basic geometric shapes will be used.
To implement the procedure from our primary aim in a test study that arises given a certain biomechanical hypothesis (see below). In order to describe the test study and the hypothesis which underlies it, we need some prerequisites and definitions, in particular, we need to define our response variable.
Our FEM (Finite Element Method) simulations were performed using Comsol Multiphysics (the successor of FEMLAB). Comsol has recently appeared in the research of the spine, e.g. for simulating a segment of the vertebral column in order to study if the stress-strain on the vertebral isthmus is affected by a bifid arc, see Quah et al. [6 (Editors) Ambrosio L., Tanner E., Biomaterials for Spinal Surgery, Woodhead Publishing; 2012]. Many other competing FEM software exist, e.g. ANSYS, used by Anburajan & Divya [7T. Antosik, "Numerical and experimental analysis of biomechanics of three lumbar vertebrae", J. Theor. Appl. Mech., vol. 37, no. 3, pp. 413-434.], to simulate lumbar spine. Other examples can be found in Aleti et al. [8D. Barnett, W. Cai, and C. Weinberger, Elasticity of Misroscopic Strustures., Lecture Notes: Stanford, .], where they use fundamental geometric figures to model the spine and among other things calculate the effect of many repetitions of certain stresses, Nédli et al. [9R.C. Brink, T.P.C. Schlösser, D. Colo, L. Vavruch, M. van Stralen, K.L. Vincken, M. Malmqvist, M.C. Kruyt, H. Tropp, and R.M. Castelein, "Anterior Spinal Overgrowth Is the Result of the Scoliotic Mechanism and Is Located in the Disc", Spine, vol. 42, no. 11, pp. 818-822.[http://dx.doi.org/10.1097/BRS.0000000000001919] [PMID: 27683977] ], where they study the compression of discs, upon applied stress to the vertebrae, and Shazly et al. [10E.N. Ebbesen, J.S. Thomsen, H. Beck-Nielsen, H.J. Nepper-Rasmussen, and L. Mosekilde, "Age- and gender-related differences in vertebral bone mass, density, and strength", J. Bone Miner. Res., vol. 14, no. 8, pp. 1394-1403.[http://dx.doi.org/10.1359/jbmr.1999.14.8.1394] [PMID: 10457272] ] & [11T. Brown, R.J. Hansen, and A.J. Yorra, "Some mechanical tests on the lumbosacral spine with particular reference to the intervertebral discs; a preliminary report", J. Bone Joint Surg. Am., vol. 39-A, no. 5, pp. 1135-1164.[http://dx.doi.org/10.2106/00004623-195739050-00014] [PMID: 13475413] ], where they model vascular blood flow and study the changes with respect to spinal chord compression. Our model construction, however, starts from DICOM stacks (underlying chest CT-scans). From this, one can create CAD-like mesh models via segmentations, and then simulate load-states and measure some appropriate output or response variable. Our chosen response variable in the present FEM study is influenced by the software chosen, i.e. Comsol. Let us now present the response variable.
First note that an applied vertical downward stress on a vertebra above the apex will render a distortion energy response, at each point on a surface that is approximately parallel to the top or bottom surface of a model of the apex. The response will, of course, look different depending upon where the stress is applied (i.e. which load-state) is used, see Definition 2.1) on the one hand and upon the curvature, wedge and tilting properties of the vertebrae on the other hand. We have chosen to perform our analysis on artificial geometric models of the apex whose interface with the intervertebral discs are modeled as planar. Our artificial geometric models have planar interface with the intervertebral discs. For our test implementation on such models, the following definitions will be used.
We define a load-state (in Comsol) to be a fixed set of specified solids in a Comsol model, together with non-void body load conditions and non-void fixed constraints. (In particular, given a load-state, the minimal conditions are fulfilled, for being able to simulate point-wise distortion energy calculations in a solid mechanics environment in Comsol).
Given a bounded planar domain, with a given coordinate grid and with pointwise calculated, non-constant values of distortion energy, for some model and a given load-state; assume that the global maximum occurs on the union of finitely many subdomains, curves and points. In the case where the global maximum is attained on a subdomain, we assign to that domain the point which constitutes the relative center of mass of the subdomain. If there is at least one subdomain, we use the point assigned to the subdomain of largest measure, with that property. That will be called the focus point. If there are several domains of the same measure the mean of their centers of mass is used. (Practical note: In cases when there is more than one subdomain and where we suspected that one of them could be the result of artefacts or inconsistencies (e.g. sharp edges) between the artificial model and the CT, manual investigation was performed and the focus point was decided ad hoc by the investigator)
Fig. (1). Consider a 3D model of a vertebral column, from an upper body CT, and let Λ be a plane intersecting the apex in the model. Assume a load-state is given. Let ℓ be the projection, on Λ, of the line (as described in Section 3.2), that describes the tilting and rotation of the hips relative to the scanner bed in the CT. The apex top response angle with respect to Λ, for the given load-state, is the angle relative to ℓ, of a line passing through (1) the relative center of mass, of the closure of the intersection of the apex with Λ, and (2) a point at which the amplitude of the distortion energy is maximal (the focus point or a point whose properties best resembles such a point, chosen ad hoc by the investigator). When the load-state, Λ and ℓ are clear from the context, we shall simply use the term apex top response angle.
In this paper, we shall only apply the term apex top response angle to the situation of an artificial geometric model, where the top face of the vertebral body in the model is planar, and , Λ in Definition 2.3, will always be assumed to contain the top face, hence we shall be working in a situation that resembles the one depicted in Fig. (1).
Now that we have defined our response variable we can formulate a hypothesis that motivates the implementation describes above as our secondary aim.
The wedge and tilting properties of the apex and nearby vertebrae, in scoliosis patients, render (for some appropriately chosen load-state) a distortion energy response, which accumulates in amplitude at an apex top response angle that varies as an analytic function of the apex vertebral rotation (for some interval).
In order for our response variable to be logical, we need to verify that there are notable variations with respect to some common load-states and patient models. In order to best suit our practical purposes we devised it to be less sensitive to shear stresses, than to transversal, almost vertical (for a standing patient) stresses. Additional shear stresses were of less interest is due to the lack of explicit surrounding tissue and ligaments in our model. To assess the latter we investigated if there was a difference between the focii of the distortion energy response when comparing load states with additional left and right shear component respectively. An additional question to address is whether the chosen response variable was affected by interchanging orthotropic material parameters with isotropic material parameters.
We calculated, for each of the basic load-states (LS), the sample correlation,
(1) |
between the vector , of apex rotations and the vector, , of the response angles for the isotropic models with the given load-state (here denotes the mean of a given vector X). We investigated whether or not a correlation was probable, and if so, whether or not that was more prominent for a particular load-state.
We point out that the main component of our work is the actual introduction of the algorithm for constructing artificial geometric models from MPR (multi-plane reconstruction) based on CT-scans. This yields a formalized method to obtain mechanically well-adapted models, which enhance and focus on certain mechanical responses which otherwise would be diminished due to small surface variations and also render shorter meshing and calculation times in FEM-analyses.
For the readers convenience we have collected in the Appendix, the details regarding distortion energy, material parameters and assignment of material properties for the orthotropic models, in practice.
All chest CT were preoperative. The scanning device was SOMATOM Definition Flash, Siemens Medical Solutions, Forchheim, Germany. The CT used 4 mm intervals. The associated software for analyzing DICOM, was Sectra, IDS7, including the application MPR.
We did this manually, based upon the method of Aaro-Dahlborn [12Q. Chen, and G. Thouas, Biomaterials: A Basic Introduction., CRC Press, .], with the requirement that the initial chosen line in the method pass the vertebral groove (for alternative methods/representatives of apex axial rotation see e.g. Lam et al. [13J.D. Currey, "The structure and mechanics of bone", J. Mater. Sci., vol. 47, pp. 41-54.[http://dx.doi.org/10.1007/s10853-011-5914-9] ]).
The method of obtaining the rotation, can be described in the following steps (we have performed the steps manually):
(1) Find an appropriate plane in 3D, that passes through an estimation of the center of mass of the apex (see Fig. 2).
(2) Find a line, in the plane from step (1), that passes through the neural groove (in the posterior of the spinal canal) such that the line divides the 2D slice of the vertebral body into two parts roughly symmetric with respect to mass.
(3) The apical rotation in the sense of Aaro-Dahlborn (with the requirement that the chosen line pass the neural groove), we shall denote it S_{AD}, is then the angle between the line in step (2), and a line passing through the same neural groove as in step (2) and also passing through the exterior mid of the sternum at the level of the plane chosen in step (1) (see Fig. 3). Though our choice of the direction of the line in step (2) was made manually, we have a posteriori, verified that our choice is compatible with the following type of symmetry about the neural groove: there is a circle in the plane chosen in step (1), such that the circle is centered at the neural groove in that plane, and defines precisely two points on the interior boundary of the sacral canal, such that the angle at the neural groove, defined by the three points has bisector, that coincides with our chosen line in step (2).
Fig. (2) The first step is to estimate an appropriate plane in 3D through the apex. |
We shall also use a representative for the tilting and rotation of sacrum (and in some sense the hip), relative to the scanner bed. Here one measures the angle, henceforth called the sacrum to table angle, between the axial horizontal with the projection of a line passing through the anterosuperiormost part of the sacro-iliac joints, in a plane which is tilted in line with the tilting of sacrum and the hips (see Fig. 4) according to the following outline:
The model extraction process from DICOM (Fig. 5), goes as follows:
(a) In the environment called MPR (Multi-Plane Reconstruction), in PACS, with the possibility of extracting tilted plane slices from the original DICOM stacks, we fix the sagittal, axial and coronal planes, such that they pass an approximate center of mass of the vertebral body, for the vertebra that we wish to model. The coordinate system will be the same in our Comsol model (i.e. based upon the scanner bed).
(b) In the fixed xz-plane we place a central line through the mid points of the upper and lower boundaries of the vertebral body (the line is depicted in the coronal plane in Fig. (5), upper left part). The height of the orthogonal projection of the vertebral body, in the slice, on this line, is used as the height of the auxiliary rotation body which we shall define in (c).
(c) Translating this line such that it passes the lower left point of the vertebra in the slice, and then drawing orthogonal lines at the lower left and upper left point of the vertebra in the slice, we obtain three sides of a preliminary rectangle. The upper and lower sides of the vertebra in the slice, yield two angles with respect to the upper and lower parts of the preliminary rectangle (the two angles are depicted in the coronal plane in Fig. (5). These angles will be used in Step 2 below for defining wedge properties.
(a) We define a trapezoid that whose one side is parallel to a vertical axis (see upper right part of Fig. 5), The height of our trapezoid will be the height defined in the last step. For the lower horizontal length, we extract a slice (in MPR) approximately parallel to the lower face of the vertebral body in 3D (depicted in upper right part of Fig. 5), and in that slice we extract the radius of an approximate circle juxtapositioned with the vertebral body in the slice. Similarly, we obtain the upper horizontal length of our trapezoid. (b) Rotating this trapezoid about its right vertical side we obtain a preliminary rotational body depicted directly below the trapezoid in Fig. (5).
Using the angles that the upper and lower sides of the vertebral body make with the preliminary rotational solid, in the xz and in the yz-planes respectively, we use Boolean operations, with auxiliary cylinders, to manifest the wedge properties in the xz and yz-planes respectively. This is depicted directly below the coronal slice in Fig. (5).
We rotate the wedged solid about the x-axis and about the y-axis respectively (it is clear that rotation about any two axes suffices to completely describe a given 3D rotation of any vector), in order to obtain the approximated tilting in 3D that we observe in the MPR (this step will affect the tensor which we need to input in Comsol, regarding material parameters, in the orthotropic case). The rotation angles are measured in the sagittal and coronal planes as depicted in Fig. (5). Fig. (6) shows a less detailed schematic of the artificial geometric vertebra construction described in the above steps.
Fig. (5) Visual aid for the described steps of the algorithm of the artificial geometric vertebra construction described stepwise in the main text. |
Fig. (6) Overview schematic of the artificial geometric vertebra construction. |
Using the approximated material parameters of Table 1, and the simplifications described below, we can summarize the flow chart of going from DICOM stacks to Comsol compatible meshed solid, in Fig. (7). FreeCAD/Meshlab), these steps are not included in the flow chart. Using these middle programs we refine, repair and convert to solid and, most importantly, convert the format to a CAD-compatible format, which is then imported into Comsol. There we attempt to mesh it. In most cases, this failed and we had to go back to the middle program for further repair and simplification. Then we also constructed a model that was simplified further using a built-in harmonic smoothing application (rendering smoother and more curved surfaces).
We investigated load-states with and without additional applied shear component (in the axial rotational sense), specifically we looked at five different load-states without axial/rotational shear, each corresponding to a reasonable situation which a spine could undergo, and then we used the same settings but with simultaneous shear (axial rotational) stress added to the processus spinosus. The latter shear addition is done both clockwise and counter-clockwise respectively. This yields a total of 5x3=15 load-states. The first four non-shear load-states where set up by fixing the lower portion of the bottom vertebra, and subdividing the upper portion of the uppermost vertebral body, into four parts (in order to do this geometrically and anatomically consistent, we defined a ’horizontal’ line to be a line parallel to the enlignment of the sacrum/hips (Section 3.2), the so called sacrum to table angle gave us this latter line. The subdivision is based upon an approximate circle determined using the ventral part of the vertebral body in the chosen slice (by ventral part we mean up to the ventral junction with the processus spinosus in the given slice). Then translate the horizontal line such that it approximately passes through the center of the above circle for the uppermost vertebra in the model (in an appropriately chosen apical slice) and we automatically obtain a vertical line in the chosen slice by using the perpendicular to our ’horizontal’. This gives us precisely four quadrants, Fig. (8). For the five non-shear load-states, the lower face of the bottommost vertebral body was used as fixed constraint, the upper face of the uppermost vertebral body was divided into four quadrants. In each of the first four load-states, precisely one of the quadrants, was subdued to downward stress. The fifth load-state is simply, fixing the lower face of the bottommost vertebra, and applying an evenly distributed stress on the upper portion of the uppermost vertebral body, i.e. all four quadrants of Fig. (8) are stressed. The additional shear analysis was done by applying shear stress on circular lateral portions, counter-clockwise and clockwise respectively (Fig. 9).
Fig. (8) Illustration of how the sacrum to table angle was used to obtain four quadrants, which were in turn used to define the load-states. |
Fig. (9) The lower face of the bottommost vertebra was fixed, and we subdivided the upper face of the uppermost vertebral body into four quadrants. In each of the first four load-states (denoted as load-states I-IV), precisely one of the quadrants I-IV, was subdued to downward stress. The fifth load-state is obtained by fixing the lower face of the bottommost vertebra, and applying an evenly distributed stress on the upper portion of the uppermost vertebral body, i.e. all four quadrants of Fig. (8) are stressed. The additional shear analysis was done by applying shear stress on circular lateral portions, counter-clockwise and clockwise respectively. The lower part of the image illustrates this for load-state V. |
Table 1 gives some basic specifications on our patient group and measurement result from the measurements made from the DICOM data using PACS, and Fig. (10) gives an illustration of the models. The similar calculations in Table 2 for the states with additional shear and for the orthotropic case where not expected to yield any better linear relations, and we did not overall see statistically significant differences in our t-tests.
Fig. (10) Illustration of the underlying models without applied stresses. |
Table 2 and Figs. (11, 12 and 13) for the apex top response angle results. Table 3 gives the results regarding the analysis of additional shear components in the loads and also the comparison between isotropic and orthotropic modeling.
Fig. (11) Apex top response angle with respect to apex rotation for the basic load states 1-5. |
Fig. (12) The load state with highest correlation coefficient |
Fig. (13) Illustration of distortion energy profiles associated to Figure 12 |
In our method investigation involving segmentation model comparison, we used a HP Elitebook 8540w, with 8Gb RAM. Meshes of type free tetrahedral coarse were used in the mesh-time comparisons analysis. In the test-runs on single apex for the manually segmented model, mesh of type coarser were used, due to a shortage of RAM in our equipment. Even so, it was clear that the calculation times were larger than for the other two methods. It makes sense to analyze the meshing time also separately (even though the simulation time of a test run often includes a meshing process) because first of all, the meshing process is completely independent of whether or not one is modeling isotropically or orthotropically with respect to the material parameters, and second of all, once a mesh is constructed, different boundary and fixed constraints will render different levels of complexity depending on how much of the rough and irregular surfaces are included for the segmentation models. In our test-runs for comparing calculation times, we did not use all three vertebrae and two discs, but rather only the apex. As fixed constraint we fixed the lower part of the vertebral body, and as boundary load we applied 1 MPa in the negative z-axis-direction, to the intersection of the top of the vertebral body, with a cylinder whose radius was adapted to the vertebral body. Fig. (14) for the set-up and execution of test-runs for the apex with isotropic material parameters. We performed the test runs with isotropic parameters and orthotropic parameters separately. The meshing times are given in Table 4 (recall that the meshing process is completely independent of whether or not one is modeling isotropically or orthotropically with respect to the material parameters). The calculation times are given in Table 5. Each of the time values in both tables was reproduced (±1 seconds) three consecutive times.
We view this as a method article, the number of patients is too small to draw statistical conclusions using correlation coefficients, but we do obtain indications on which load-states could be of interest for further experiment set-ups. Also we do not necessarily expect to find linear relations between our response variable and apex rotation but an analytic relation could exist.
The main pros of our method involving artificial geometric models is of course that it enhances the semi-global mechanical responses to applied stresses and that it requires less calculation times due largely to less complicated interfaces. The main cons are that local details are lost in the process (which for example should be present in the mechanical analysis of screws and rods interacting with the vertebrae or in the study of morphology). Note that the Cobb angle is a parameter that is not local to the apex, indeed the neighbours of the apex and their relation to it, do not necessarily reflect the severity of the total coronal curvature. On the other hand, the apex rotation, as we have used in this paper, is strictly local and confined to the apex itself, and thereby a more logical choice of deviation parameter (for our purposes). Here are some restrictions that are made in the pre-processing and processing stages:
For mechanical analysis of manual segmentation models, we would need to adjust or replace our response variable, because the local irregularities should give false peaks if one only looks at slices of the original 3D-data. It would then possibly be more appropriate to extract some thin volume which one subdivides and performs integration of distortion energy over the elements obtained by the subdivision.
We expected, due to our choice of variable, that additional shear stresses would not render a linear relation between our response variable and apex rotation, as opposed to strictly transversal (rather vertical) stresses. This was verified by pairwise t-tests. The reason additional shears were of less interest is due to the lack of explicit surrounding tissue and ligaments in our model. However, we had no prior expectation on whether or not there would be a statistically significant difference between isotropic and orthotropic modeling (which it turned out in general not to be the case). Obviously, the small amount of patient models means that we cannot draw any conclusion about the linearity between our response variable and apex rotation, but we can at best say that we have obtained indications that a correlation could exist between apex rotation and the response variable, for some appropriately chosen load-state, and our results indicate what requisites such a load-state should have.
We note that all our patients had right-convex scoliosis and a slight local lordosis at the apex. The load-state rendering the best correlation coefficient is load-state 2, and it is the only among the 5 basic load-states with the following properties:
(1) In the coronal plane, the stress counteracts the curvature (in the sense that it lies on the thicker part of the wedged slice of the apex).
(2) In the sagittal plane, the stress counteracts the curvature (in all cases a slight lordosis). If this observation can be reinforced in the future with more substantial evidence, we do not have a good explanation for part (2) of the observation.
We introduce an algorithm for constructing artificial geometric models from CT-scans and thereby a formalized way to obtain models, which enhance and focus on certain mechanical responses which otherwise would be diminished due to small surface variations. Calculation and meshing times appear much shorter. We chose a particular response variable based upon distortion energy distribution, which our pairwise t-tests show is non-trivial, and is foremost sensitive to the type of transversal stresses we are interested in. The reason additional shears were of less interest is due to the lack of explicit surrounding tissue and ligaments in our model. The method is inappropriate for analyses where the local roughness and irregularities of surfaces are wanted features. We view this as a method article, the number of patients is too small to draw statistical conclusions using correlation coefficients, but we do obtain indications of e.g. which load states are promising candidates with regards to the most common patient specifications. The load-state for which we obtained highest correlation between apex rotation and our chosen response angle, was one where the stress, locally, counteracted the coronal curvature and simultaneously counteracted the sagittal curvature. In general, our model does not render different results for isotropic and orthotropic material modeling respectively. By creating three different types of models for a single patient we illustrate the gain in mechanical simplicity, but a loss in detail using our artificial models. We have pointed out in the article that, if we would like to perform analogous analysis on manually segmentation models in future work, an integration process with appropriately chosen volume, at the interesting faces of the apex, and subdivision of that volume, could be one possible way to substitute the apex top response angle.
The study was funded by the Swedish government in terms of ALF-grant (Event No:LIO-608221).
The study was approved by the regional ethics committee in Linköping, Sweden (2012/366-31).
Animals did not participate in this research. All human research procedures followed were in accordance with the ethical standards of the committee responsible for human experimentation (institutional and national), and with the Helsinki Declaration of 1975, as revised in 2008.
Consent for publication is obtained.
The authors declare no conflict of interest, financial or otherwise.
Declared none.
The failure theory of von Mises, based upon maximum distortion energy, roughly assumes that some notion of ’failure’, of a medium under stress, occurs when the maximum distortion energy at a reference point, in a stressed material reaches the value of the maximum distortion energy at the failure point in simple uniaxial tension (this failure stress is denoted by σ_{failure tensile}), and the distortion energy is given by the left hand side of the following equation,
(A1) |
Where σ_{i}, i = 1, 2, 3, denote the principal stresses. Akin [14V.H. Frankel, and M. Nordin, Basic Biomechanics of the Musculoskeletal System, 3 ed., 2001. ], p.267, defines the von Mises stress as the square root of, one half of the left hand side of equation 2, namely,
(A2) |
and relates the von Mises stress to the distortion energy. Although the von Mises energy definition was constructed based on the isotropic case, one can use the same definition for any case when the principle stresses are calculated, and it is then an approximative entity for describing distortion energy matters.
For the Comsol model one must have certain physical constants and parameters (or approximations of them) at hand. Bones and other human tissue are very complex materials chemically and physically. Bone is in general a connective tissue consisting of cells embedded in a mineralized matrix including collagen fibers which are impregnated with calcium phosphate mainly as hydroxyapatite, as well as carbonate, citrate, sodium, and magnesium. Bone is in general composed of about 75% inorganic material and 25% organic material. Different types of bones vary in physical properties, e.g. the diaphyses of the tibia and femur have a rather large marrow, and a complicated Haverian system of canals, whereas such features are much less prominent in e.g. the frontal cranium. One of the features of certain bones, such as diaphyseal bones and also to smaller extent, the vertebral bodies that encourages an orthotropic modeling, is the fact that osteones have concentric lamellae that induce an intrinsic main/longitudinal axis. The densities of bone and vertebral discs are sometimes required inputs for defining a material in FEM models, and in such case we have used 1908 kg/m3 (the software built-in value of density of generalized bone). For the intervertebral discs, we use, following Chen & Thouas [15J.M. Mac-Thiong, F.M. Pinel-Giroux, J.A. de Guise, and H. Labelle, "Comparison between constrained and non-constrained Cobb techniques for the assessment of thoracic kyphosis and lumbar lordosis", Eur. Spine J., vol. 16, no. 9, pp. 1325-1331.[http://dx.doi.org/10.1007/s00586-007-0314-1] [PMID: 17426991] ], p.400, the density of dry collagen to be approximately 1300 kg/m3, and we approximate the density equivalent composition as 60% water and 40% collagen, (in fact, both the composition for the annulus fibrosus and nucleus pulposus contain also proteoglycan and third main component, thus we have made an approximation of using the dry weight of collagen to approximate that of proteoglycan). Even after approximations regarding mechanical and physical properties there remains other complications such as age and gender, see Ebbesen et al. [16F. Hefti, "Pathogenesis and biomechanics of adolescent idiopathic scoliosis (AIS)", J. Child. Orthop., vol. 7, no. 1, pp. 17-24.[http://dx.doi.org/10.1007/s11832-012-0460-9] [PMID: 24432054] ]. We mention some research articles including material constant estimation, see e.g. Antosik et al. [17F. Irgens, Continuum Mechanics., Springer, .], Kurutz et al. [18T.M. Keaveny, Standard handbook of biomedical engineering and design., McGraw-Hill, .] and Nédli et al. [19C. Kittel, Introduction to solid state physics., 7th edJohn Wiley & Sons, .7th ed]. Other sources in which there appear material constant estimates are Keaveny et al. [20M. Kurutz, "Age-sensitivity of time-related in vivo deformability of human lumbar motion segments and discs in pure centric tension", J. Biomech., vol. 39, no. 1, pp. 147-157.[http://dx.doi.org/10.1016/j.jbiomech.2004.10.034] [PMID: 16271599] ], Ch.8, Currey [21M. Kurutz, Finite Element Modeling of the Human Lumbar Spine, Finite Element Analysis (David Moratal (Ed), ISBN 978-953-307-123-7), 2010 ] and Frankel & Nordin [22G.C. Lam, D.L. Hill, L.H. Le, J.V. Raso, and E.H. Lou, "Vertebral rotation measurement: a summary and comparison of common radiographic and CT methods", Scoliosis, vol. 3, no. 16, p. 16.[http://dx.doi.org/10.1186/1748-7161-3-16] [PMID: 18976498] ]. Fracture stress of a general bone has in some contexts been estimated to be 14±25% (MPa), see Vable [23P. Nédli, Contact and No-Compression Analysis of a Human Spine Segment: Theory, Method and Parametric Investigation., COMSOL Users Conference Grenoble, .], p.3.
However it is commonly known that bones behave anisotropically and under tension, and viscoelastically (i.e. the stress-strain curves depend upon the rate of load application), see e.g. Frankel & Nordin [22G.C. Lam, D.L. Hill, L.H. Le, J.V. Raso, and E.H. Lou, "Vertebral rotation measurement: a summary and comparison of common radiographic and CT methods", Scoliosis, vol. 3, no. 16, p. 16.[http://dx.doi.org/10.1186/1748-7161-3-16] [PMID: 18976498] ]. The primary sources however, used in the collecting of material constants and insight to the material properties for us have been: White [25C. Quah, M.S. Yeoman, A. Cizinauskas, K.C. Cooper, N.S. Peirce, D.S. McNally, and B.M. Boszczyk, "Finite element investigation of the effect of a bifid arch on loading of the vertebral isthmus", Spine J., vol. 14, no. 4, pp. 675-682.[http://dx.doi.org/10.1016/j.spinee.2013.08.040] [PMID: 24268389] ], Kurutz [18T.M. Keaveny, Standard handbook of biomedical engineering and design., McGraw-Hill, .], Wilke et al. [26H. Schmidt, A. Kettler, F. Heuer, U. Simon, L. Claes, and H.J. Wilke, "Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading", Spine, vol. 32, no. 7, pp. 748-755.[http://dx.doi.org/10.1097/01.brs.0000259059.90430.c2] [PMID: 17414908] ], Stemper et al. [27T. Shazly, "Effect of Spinal Cord Compression on Local Vascular Blood Flow and Perfusion Capacity", Proceedings of the 2012 Comsol Conference, .], Schmidt et al. [24P.O. Newton, Idiopathic Scoliosis: The Harms Study Group Treatment Guide., Thieme, .], Ambrosio & Tanner [28T. Shazly, Assessment of Anterior Spinal Artery Blood Flow following Spinal Chord Injury.PLoS ONE , vol. 9, no. 9, pp. 1-7. e108820, doi:10.1371/journal.pone.0108820] (and indirectly Brown, Hansen and Yorra [29B.D. Stemper, D. Board, N. Yoganandan, and C.E. Wolfla, "Biomechanical properties of human thoracic spine disc segments", J. Craniovertebr. Junction Spine, vol. 1, no. 1, pp. 18-22.[http://dx.doi.org/10.4103/0974-8237.65477] [PMID: 20890410] ]), and the references found in these later articles. Interpretation of such data requires a prescribed main axis (in the case of our cited articles, the z-axis) in the coordinate system for which the data holds true, and a basic knowledge of the anatomy and structure of bone makes it clear that the main axis can, for the epiphyseal case, be chosen along the Haverian longitudinal axis, and for vertebra, as one passing through the center of the vertebral body and approximately tangential to the spinal cord. We have chosen to mainly use the constants gathered from different literature in Schmidt et al. [24P.O. Newton, Idiopathic Scoliosis: The Harms Study Group Treatment Guide., Thieme, .]. Note that such approximations impose the difficulty that osteoporosis, if not taken into account, may distort the authenticity of the simulation.
The feature we use in Comsol, starts with applied force, applied loads and fixed constraints. The output is a version of distortion energy at each point (we use stationary analysis, as opposed to time-dependent). The calculations performed by Comsol thus requires the input of material parameters in terms of either compliance or elasticity components. Assuming linearly elasticity, there exists a tensor (see e.g. Irgens [30I.A.F. Stokes, Mechanical effects on skeletal growth, J Musculoskel Neuron Interact, 2(3), (2002), 277-31. Vable M., Mechanics of Materials., Michigan Technological University, .]), S, let us denote its components by s_{ijkl,} interpreted as a fourth order tensor, which acts on strain tensors (ε, which are second order tensors in their own right) and yields the associated stress tensor (σ, also a tensor of second order) and is sometimes referred to as the elasticity tensor. The elasticity tensor is the inverse (in the appropriate sense) of the so called compliance tensor, let us denote it C. Both the strain and stress tensors are symmetric, thus each is completely determined by 6 of its components (namely σ_{ij} = σ_{ji} for i≠j, so we only need to know the upper triangle part of a matrix representation of σ), and we can represent C by 36 components in matrix form, using the so-called Voigt map of indices:
(A3) |
(A4) |
(A5) |
which induces,
(A6) |
The tensors C (and analogously S) can be shown (see e.g. Kittel [31M. Vable, Mechanics of Materials, Michigan Technological University, 2009] or Zehnder [32M.M. Panjabi, and A.A. White III, "Basic biomechanics of the spine", Neurosurgery, vol. 7, no. 1, pp. 76-93.[http://dx.doi.org/10.1227/00006123-198007000-00014] [PMID: 7413053] ], p.15) to be symmetric in the sense that,
(A7) |
The relation between the arrays S and C can be given in tensor form, is given by (see e.g. Barnett et al. [33H.J. Wilke, S. Krischak, and L.E. Claes, "Formalin fixation strongly influences biomechanical properties of the spine", J. Biomech., vol. 29, no. 12, pp. 1629-1631.[http://dx.doi.org/10.1016/S0021-9290(96)80016-9] [PMID: 8945663] ], p.15),
(A8) |
where δ_{ij} for i=j, and zero otherwise. The right hand side of the last equation has a 6x6 representation, when applying the Voigt map to the index pairs (i,j) and (n,m).
Let η_{1},η_{2},η_{3}, be indices obtained by applied the Voigt map to the index pairs (i,j), (k,l) and (n,m) respectively, and C_{η1η2}, S_{η1η2} let denote the tensor components after having transformed the indices. Then the left hand side of equation (A8) can (using the symmetry relations of equation (A7)) be described as,
(A9) |
Combining the last two representations we obtain 36 equations which can be represented in matrix form according to (note that some text books prefer to use c for the elements of the elasticity tensor, we have chosen to use c for the elements of the compliance tensor),
Which is equivalent to,
In other words, as square matrices C^{voigt}, S^{voigt}, satisfy, (C^{voigt})^{-1} = S^{voigt} For an orthotropic material, we can (given the Young/shear moduli and Poisson ratios: Y_{xx}, Y_{yy}, Y_{zz}, Y_{yz}, Y_{xz}, Y_{xy}, V_{xy}, V_{yz}, V_{xz} for appropriate choice of basis, write the compliance matrix, in Voigt notation as,
which will be symmetric. In our case, based upon the references we used, we specifically consider the case Y_{xx} = Y_{yy}, Y_{yz}, Y_{xz} (such a model is sometimes called transversely orthotropic) and also V_{yz} = V_{xz}. Note that, by symmetry, V_{xy} = V_{yx}, V_{zy} = V_{zx} and (again in the special case within which we work) If C = [C_{ijkl}] is the compliance array (34 elements), we can for a given rotation in 3D, given by say the matrix R; use the relation (see e.g. Irgens [29B.D. Stemper, D. Board, N. Yoganandan, and C.E. Wolfla, "Biomechanical properties of human thoracic spine disc segments", J. Craniovertebr. Junction Spine, vol. 1, no. 1, pp. 18-22.[http://dx.doi.org/10.4103/0974-8237.65477] [PMID: 20890410] ], p.88),
(A10) |
where can be interpreted as the compliance matrix of an orthotropic material with respect to the new basis R(e_{1}, e_{2}, e_{3})^{T} We can then use Voigt notation to extract the compliance and elasticity array for a ’tilted’ orthotropic material using Comsol. As an example, rotation about the x-axis of Ø_{1} degrees and rotation about the y-axis of Ø_{2} degrees respectively are given (if we choose a setting where x increases to the right, y points out of the screen and z increases upward, thereby the following matrices will for positive angles give clockwise rotation, when the axis of rotation points toward observer) by the orthogonal transformation matrices,
The combined transformation matrix becomes, R = R^{x, Ø1}R^{y, Ø2} Using the above formulae, a Matlab code was written in order to calculate the input elasticity matrix elements for tilted vertebrae in the orthotropic models.
[1] | S. Aaro, and M. Dahlborn, "Estimation of vertebral rotation and the spinal and rib cage deformity in scoliosis by computer tomography", Spine, vol. 6, no. 5, pp. 460-467.[http://dx.doi.org/10.1097/00007632-198109000-00007] [PMID: 7302680] |
[2] | American Association of Neurological Surgeons, http:// www.aans.org/ Patient%20Information/ Conditions%20and%20Treatments/ Scoliosis.aspx |
[3] | M. Anburajan, and V. Divya, "Finite Element Analysis of Human Lumbar Spine", 2011 International Conference on Electronics Computer Technology, . |
[4] | J.E. Akin, Finite Elements for Analysis and Design., Academic Press, . |
[5] | H. Aleti, and I.A. Motaleb, "Study of Stress and Strain in Human Spine", IEEE International Conference, pp. 184-9.[http://dx.doi.org/10.1109/EIT.2014.6871759] |
[6] | (Editors) Ambrosio L., Tanner E., Biomaterials for Spinal Surgery, Woodhead Publishing; 2012 |
[7] | T. Antosik, "Numerical and experimental analysis of biomechanics of three lumbar vertebrae", J. Theor. Appl. Mech., vol. 37, no. 3, pp. 413-434. |
[8] | D. Barnett, W. Cai, and C. Weinberger, Elasticity of Misroscopic Strustures., Lecture Notes: Stanford, . |
[9] | R.C. Brink, T.P.C. Schlösser, D. Colo, L. Vavruch, M. van Stralen, K.L. Vincken, M. Malmqvist, M.C. Kruyt, H. Tropp, and R.M. Castelein, "Anterior Spinal Overgrowth Is the Result of the Scoliotic Mechanism and Is Located in the Disc", Spine, vol. 42, no. 11, pp. 818-822.[http://dx.doi.org/10.1097/BRS.0000000000001919] [PMID: 27683977] |
[10] | E.N. Ebbesen, J.S. Thomsen, H. Beck-Nielsen, H.J. Nepper-Rasmussen, and L. Mosekilde, "Age- and gender-related differences in vertebral bone mass, density, and strength", J. Bone Miner. Res., vol. 14, no. 8, pp. 1394-1403.[http://dx.doi.org/10.1359/jbmr.1999.14.8.1394] [PMID: 10457272] |
[11] | T. Brown, R.J. Hansen, and A.J. Yorra, "Some mechanical tests on the lumbosacral spine with particular reference to the intervertebral discs; a preliminary report", J. Bone Joint Surg. Am., vol. 39-A, no. 5, pp. 1135-1164.[http://dx.doi.org/10.2106/00004623-195739050-00014] [PMID: 13475413] |
[12] | Q. Chen, and G. Thouas, Biomaterials: A Basic Introduction., CRC Press, . |
[13] | J.D. Currey, "The structure and mechanics of bone", J. Mater. Sci., vol. 47, pp. 41-54.[http://dx.doi.org/10.1007/s10853-011-5914-9] |
[14] | V.H. Frankel, and M. Nordin, Basic Biomechanics of the Musculoskeletal System, 3 ed., 2001. |
[15] | J.M. Mac-Thiong, F.M. Pinel-Giroux, J.A. de Guise, and H. Labelle, "Comparison between constrained and non-constrained Cobb techniques for the assessment of thoracic kyphosis and lumbar lordosis", Eur. Spine J., vol. 16, no. 9, pp. 1325-1331.[http://dx.doi.org/10.1007/s00586-007-0314-1] [PMID: 17426991] |
[16] | F. Hefti, "Pathogenesis and biomechanics of adolescent idiopathic scoliosis (AIS)", J. Child. Orthop., vol. 7, no. 1, pp. 17-24.[http://dx.doi.org/10.1007/s11832-012-0460-9] [PMID: 24432054] |
[17] | F. Irgens, Continuum Mechanics., Springer, . |
[18] | T.M. Keaveny, Standard handbook of biomedical engineering and design., McGraw-Hill, . |
[19] | C. Kittel, Introduction to solid state physics., 7th edJohn Wiley & Sons, .7th ed |
[20] | M. Kurutz, "Age-sensitivity of time-related in vivo deformability of human lumbar motion segments and discs in pure centric tension", J. Biomech., vol. 39, no. 1, pp. 147-157.[http://dx.doi.org/10.1016/j.jbiomech.2004.10.034] [PMID: 16271599] |
[21] | M. Kurutz, Finite Element Modeling of the Human Lumbar Spine, Finite Element Analysis (David Moratal (Ed), ISBN 978-953-307-123-7), 2010 |
[22] | G.C. Lam, D.L. Hill, L.H. Le, J.V. Raso, and E.H. Lou, "Vertebral rotation measurement: a summary and comparison of common radiographic and CT methods", Scoliosis, vol. 3, no. 16, p. 16.[http://dx.doi.org/10.1186/1748-7161-3-16] [PMID: 18976498] |
[23] | P. Nédli, Contact and No-Compression Analysis of a Human Spine Segment: Theory, Method and Parametric Investigation., COMSOL Users Conference Grenoble, . |
[24] | P.O. Newton, Idiopathic Scoliosis: The Harms Study Group Treatment Guide., Thieme, . |
[25] | C. Quah, M.S. Yeoman, A. Cizinauskas, K.C. Cooper, N.S. Peirce, D.S. McNally, and B.M. Boszczyk, "Finite element investigation of the effect of a bifid arch on loading of the vertebral isthmus", Spine J., vol. 14, no. 4, pp. 675-682.[http://dx.doi.org/10.1016/j.spinee.2013.08.040] [PMID: 24268389] |
[26] | H. Schmidt, A. Kettler, F. Heuer, U. Simon, L. Claes, and H.J. Wilke, "Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading", Spine, vol. 32, no. 7, pp. 748-755.[http://dx.doi.org/10.1097/01.brs.0000259059.90430.c2] [PMID: 17414908] |
[27] | T. Shazly, "Effect of Spinal Cord Compression on Local Vascular Blood Flow and Perfusion Capacity", Proceedings of the 2012 Comsol Conference, . |
[28] | T. Shazly, Assessment of Anterior Spinal Artery Blood Flow following Spinal Chord Injury.PLoS ONE , vol. 9, no. 9, pp. 1-7. e108820, doi:10.1371/journal.pone.0108820 |
[29] | B.D. Stemper, D. Board, N. Yoganandan, and C.E. Wolfla, "Biomechanical properties of human thoracic spine disc segments", J. Craniovertebr. Junction Spine, vol. 1, no. 1, pp. 18-22.[http://dx.doi.org/10.4103/0974-8237.65477] [PMID: 20890410] |
[30] | I.A.F. Stokes, Mechanical effects on skeletal growth, J Musculoskel Neuron Interact, 2(3), (2002), 277-31. Vable M., Mechanics of Materials., Michigan Technological University, . |
[31] | M. Vable, Mechanics of Materials, Michigan Technological University, 2009 |
[32] | M.M. Panjabi, and A.A. White III, "Basic biomechanics of the spine", Neurosurgery, vol. 7, no. 1, pp. 76-93.[http://dx.doi.org/10.1227/00006123-198007000-00014] [PMID: 7413053] |
[33] | H.J. Wilke, S. Krischak, and L.E. Claes, "Formalin fixation strongly influences biomechanical properties of the spine", J. Biomech., vol. 29, no. 12, pp. 1629-1631.[http://dx.doi.org/10.1016/S0021-9290(96)80016-9] [PMID: 8945663] |
[34] | M. Zehnder, Atomistic Simulation of the Elasticity of Polymers. Diss. ETH No 12497, 1997 |