System and Method for Simulating Intervertebral Discs

Abstract

The present invention provides a method of simulating an intervertebral disc, comprising: automatically segmenting a medical image of a spine based on a pre-trained model so as to obtain different intervertebral disc segmentations; and biomechanical modelling under multiple motions for the different intervertebral disc segmentations so as to obtain dynamic properties of the intervertebral disc.

Description

System and Method for Simulating Intervertebral DiscsTECHNICAL FIELD

[0001] This invention generally relates to the field of medicine, specifically relates to a system and method for simulating intervertebral discs (IVD) , more specifically, relates to a system and method for simulating intervertebral discs (IVD) for optimizing material properties (incompressive, anisotropic, heterogeneous, etc) of a point clouds-based simulation to match those of IVD.BACKGROUND

[0002] With the development of medicine, body implants have been widely used clinically. At present, both clinicians and implant designers lack tools for generating highly accurate and bio-fidelity predictions of implant performance in bone tissue, leading to a higher likelihood of complications in the case of clinical use, and impeding the progress of the design of safer and more effective implants.

[0003] In recent decades, diagnostic imaging in orthopaedics has greatly improved in terms of its capacity to deliver high fidelity 3D images of a patient’s muscoloskeletal structure and anatomy, but these images are generally then interpreted in the traditional manual manner by surgeons and radiologists. While the skill of a well-trained clinician is irreplaceable in terms of its role in diagnosis and treatment of intervertebral disc diseases, using computational methods such as mechanical simulation can greatly enhance the decision-making process of the clinician by providing them with reference information.

[0004] The existing simulation systems fall into two categories, specialized finite-element simulation systems for body tissues simulation, which do not accurately predict large deformation, damage evolution, and cracking propagation of IVD due to physiological loading, and general-purpose  discrete-element simulation systems that have not been developed into systems capable of representing soft materials such as IVD. Existing software packages cannot achieve this result due to their reliance on inappropriate modelling and simulation methods and / or lack of a built-in implant performance prediction system. Particularly, current existing simulation method cannot represent the material properties of the soft tissues, due to the complex structure and material properties of the soft tissues.SUMMARY

[0005] The following presents a simplified summary of the invention in order to provide a basic understanding of some aspects of the invention. This summary is not an extensive overview of the invention. It is intended to neither identify key or critical elements of the invention nor delineate the scope of the invention. Rather, the sole purpose of this summary is to present some concepts of the invention in a simplified form as a prelude to the more detailed description that is presented hereinafter.

[0006] In one aspect, the present invention provides a method of simulating an intervertebral disc, comprising:

[0007] automatically segmenting a medical image of a spine based on a pre-trained model so as to obtain different intervertebral disc segmentations; and

[0008] biomechanical modelling under multiple motions for the different intervertebral disc segmentations so as to obtain dynamic properties of the intervertebral disc.

[0009] Preferably, the motions comprise a force loading.

[0010] Preferably, the force loading is simulated with a meshfree method.

[0011] Preferably, the meshfree method comprises point clouds-based simulation.

[0012] Preferably, the intervertebral disc segmentations comprise a normal disc, a bulging degenerative disc, a herniated disc and / or a vertebra.

[0013] Preferably, the intervertebral disc comprises an annulus fibrosis and a nucleus pulposus.

[0014] Preferably, the nucleus pulposus is simulated using a Neo-Hookean model.

[0015] Preferably, a Neo-Hookean strain energy potential of the Neo-Hookean model is:

[0016] where U is the Neo-Hookean strain energy potential; C10 is material constants characterizing the deviatoric deformation of the material; and D1 is material incompressibility parameter;  is the first deviatoric strain invariant; Jel is the elastic volume ratio.

[0017] Preferably, the annulus fibrosis is simulated using a Holzapfel-Gasser-Ogden model to represent the nonlinear mechanical behavior of collagen fibers.

[0018] Preferably, a Holzapfel-Gasser-Ogden strain energy potential of the Holzapfel-Gasser-Ogden model is as follows:

[0019] wherein,  and

[0020] U is the Holzapfel-Gasser-Ogden strain energy potential; C10 is material constants characterizing the deviatoric deformation of the material; D1 is material incompressibility parameter; k1 is a dimensionless parameter; k2 is a stress-like parameter to be determined from mechanical tests of the tissue; κ is structure parameters to be determined from histological data of the tissue; N is the number of families of fibers;  is the first deviatoric strain invariant; Jel is the elastic volume ratio.

[0021] Preferably, when biomechanical modelling under multiple motions for the different intervertebral disc segmentations, considering a fiber mechanics model expressed as follows:

[0022] wherein, σf and εf correspondingly represent fiber stress and strain,  representing elasticity and creep, df is a damage variable, 0≤df≤1, 0 representing no damage to the fiber, and 1 representing completely tear.

[0023] Preferably, the pre-trained model is a deep learning model, which are trained by using the samples of medical images of spines.

[0024] Preferably, the model is trained by the following steps:

[0025] Collecting the medical images of spines;

[0026] Segmenting, masking and labeling intervertebral disc, nucleus, and spinal canal of the images, to obtain the samples of medical images of spines;

[0027] training the deep learning model by using the samples of medical images of spines.

[0028] In another aspect, the present invention provides a computer-readable storage medium comprising a program, the program is executed to perform the method of simulating an intervertebral disc of the present invention.

[0029] In yet another aspect, the present invention provides a system for simulating an intervertebral disc, comprising a processor and a storage, wherein, the processor is configured to perform the method of simulating an intervertebral disc by execute the executable commands in the storage.

[0030] To the accomplishment of the foregoing and related ends, the invention comprises the features hereinafter fully described and particularly defined in the claims. The following description and the annexed drawings set forth in detail certain illustrative aspects and implementations of the invention. These are indicative, however, of but a few of the various ways in which the principles of the invention may be employed. Other objects, advantages and novel features of the invention will become apparent from the following detailed description of the invention when considered in conjunction with the drawings.

[0031] BRIEF SUMMARY OF THE DRAWINGS

[0032] Fig. 1 is a flowchart of simulating an intervertebral disc according to the  embodiment of the present invention.

[0033] Figs. 2A and 2B show the preliminary results of the intervertebral disc morphology detection.

[0034] Figs. 3A-3C show an experimental design of IVD force loading.

[0035] Fig. 4 shows an example of the compression experimental chamber for assessing disc dynamic viscoelastic properties.

[0036] Fig. 5 shows a flowchart of the research plan including experiments, computer modelling and clinical validation.DETAILED DESCRIPTION

[0037] Disclosed herein is a method to simulate spine soft tissue material properties under loading, wherein, the material properties of a point clouds-based simulation is optimized to match those of spine. Please refer to a flowchart as shown in Fig. 1, the method of simulating an intervertebral disc includes the following steps.

[0038] Step S1: Automatically segmenting various spinal structure utilizing a conventional supervised method or a unique unsupervised method with a voting mechanisms for accurate segmentation of different anatomical structures for further generation of the accurate 3D model of the spine with vertebrae, IVD and nerves.

[0039] As to the unique unsupervised method, please refer to US Non-Provisional Application Serial No. 17 / 854, 230 claiming priority of US application No. 63 / 217, 157, the entirety of which is incorporated herein by reference.

[0040] In order to generate the appropriate geometry of the IVD and to automatically assess the spine alignment, a few substeps are presented as follows:

[0041] i) Data collection and training data labelling.

[0042] The data can be collected through using medical imaging of spines under ethical approval. For data labelling task, three or more spine surgeons and radiologists can be enrolled for manually masking objects. The intervertebral  disc, nucleus, and spinal canal would be separately segmented and masked. The masks can be stored as different labels for automatic segmentation model training.

[0043] ii) Implementation of deep learning model for intervertebral disc and nucleus segmentation.

[0044] For automatically intervertebral disc recognition, the methods disclosed in the filed US patent application (US NonProvisional Application Serial No. 17 / 854,230 claiming priority of US application No. 63 / 217,157) is incorporated herein by reference. Figs. 2A and 2B show the preliminary results of the intervertebral disc morphology detection, wherein, Fig. 2A shows the original image of the spine, and Fig. 2B shows the image of the spine being segmented, wherein the normal disc, the bulging degenerative disc, the herniated disc and the vertebra are clearly marked.

[0045] Step S2: Biomechanical modelling under multiple motions for different segmentations obtained from Step 1, such as for different disc morphologies of Fig. 2B, to develop the specific algorisms (such as novel iterative methods) and determine specific material property parameters of annulus fibrosis (AF) and nucleus pulposus (NP) for simulating disc herniation and disease progression at different disease stage in response to different biomechanical influence.

[0046] As well known for a person skilled in the art, IVD has two major parts, i.e., annulus fibrosis which is fiber and nucleus which comprises mainly water and is in-compressive.

[0047] For numerical modelling, the nucleus of the IVD is simulated using the Neo-Hookean model to represent the ability of large deformations and incompressibility. The form of the Neo-Hookean strain energy potential is expressed in Equation (1) :

[0048] where U is the strain energy per unit of reference volume (i.e., strain energy potential) ; C10 is material constants characterizing the deviatoric deformation  of the material; and D1 is material incompressibility parameter;  is the first deviatoric strain invariant; Jel is the elastic volume ratio.

[0049] In terms of AF, the substance is reinforced by collagen fibers with alternating orientations. The distribution of collagen fibers generates the material anisotropy and heterogeneity. Therefore, the AF material is simulated using the Holzapfel-Gasser-Ogden model to represent the nonlinear mechanical behavior of collagen fibers. The Holzapfel-Gasser-Ogden (HGO) model is one of the most popular anisotropic hyperelastic models. Generally, those skilled may modify or improve the model by adjusting the parameters according to some initial experimental results. The form of the Holzapfel-Gasser-Ogden strain energy potential refers to the research of Holzapfel, Gasser, and Ogden (2000) and Gasser, Ogden, and Holzapfel (2006) for modeling arterial layers with distributed collagen fiber orientations:

[0050] with where U is the strain energy per unit of reference volume (i.e., strain energy potential) ; C10 is material constants characterizing the deviatoric deformation of the material; D1 is material incompressibility parameter; k1 is a dimensionless parameter; k2 is a stress-like parameter to be determined from mechanical tests of the tissue; and κ is structure parameters to be determined from histological data of the tissue; N is the number of families of fibers;  is the first deviatoric strain invariant; Jel is the elastic volume ratio. The first two terms in the Equation (2) represent the distortional and volumetric contributions of the AF ground material without fibers, and the third term represents the contributions from the collagen fibers with different orientations. And a basic assumption is that the fibers can only withstand tension, and would buckle under compression. If the mechanical behavior of AF material is too complex to be fitted by adjusting parameters, those skilled may modify or add a correction term in expression (2) . In fact, by comparing  Equations (1) and (2) , it can be seen that the HGO model itself is a modification of the Neo-Hookean model. The original HGO model was developed for anisotropic arteries, where the material is represented using 2 families of collagen fibers embedded into an isotropic matrix.

[0051] Therefore, the inventor improves the Holzapfel-Gasser-Ogden model by adjusting of parameters according to some initial experimental results, to develop a more suitable material constitutive relationship for AF of the IVD to simulate the nonlinear mechanical behavior of collagen fibers. In addition, AF of the IVD is complex that long-term creep and tear need to be considered. Therefore, the fiber mechanics model to be considered in the study is expressed in Equation (3) .

[0052] Wherein, U is the HGO strain energy potential in Equation (2) , σf and εf correspondingly represent fiber stress and strain,  representing elasticity and creep, indicates that fiber stress and strain is associated with strain rate. df is the damage variable, 0≤df≤1, 0 representing no damage to the fiber, and 1 representing completely tear. This equation (3) is the expression of stress-strain relation.

[0053] Due to the large deformation of the IVD under force loading, it is proposed to use meshfree method such as point clouds to simulate the force loading and deformation, due to the adaptability to large deformations of Point Clouds. In the field of numerical analysis, meshfree methods are those do not require connection between nodes of the simulation domain, i.e. a mesh, but are rather based on interaction of each node with all its neighbors. As a consequence, original extensive properties such as mass or kinetic energy are no longer assigned to mesh elements but rather to the single nodes. Point clouds method is also one of meshfree methods which share a common feature: the approximation of unknowns in the PDE is constructed based on scattered points without mesh connectivity.

[0054] Biomechanical testing

[0055] In order to verify the effectiveness of the above simulation, the inventor conducts a biomechanical testing to examine the association between different morphology of degenerative intervertebral disc (such as disc height deduction and nucleus migration mechanisms, etc. ) and altered biomechanical characteristics using material testing system.

[0056] material testing and Digital image correlation

[0057] Selections of the typical stress-deformation mode as shown in Figs. 3A-3C are made based on human spine biomechanics. Depending on the position of the gravity line / axis of gravity, the loading vector through the IVD may be variable. Figs. 3A-3C show an experimental design of IVD force loading to mimic the physiological and pathological conditions of the IVD, wherein, F represents longitudinal force, P represents horizontal pressure, e represents eccentric distance, Δ represents absolute displacement of the compressive plate to the IVD in the axial direction, θ represents angle of the superior compressive plate to the IVD in comparison with the initial position , Δ1 and Δ2 respectively represent longitudical and horizontal absolute displacements. In axial compression (please see Fig. 3A) , the load is through the center of the disc which occurs in a flatback when the gravity line falls through the center of the spinal segments. In such instance, resultant hoop stresses should be equal circumferentially. In eccentric compression (please see Fig. 3B) , the load may fall posteriorly in a lordotic segment creating asymmetric loading and increased posterior herniations. Similarly, in a kyphotic segment, the load may fall anteriorly and lead to anterior herniations. Compressive shearing may occur in cases of spondylolisthesis where the spinal segment is unstable causing a shearing moment in addition to axial loads (please see Fig. 3C) .

[0058] The axial compression of the IVD and eccentric compression can be tested using the material testing system. Fig. 4 shows an example of the  compression experimental chamber for assessing disc dynamic viscoelastic properties. In this setup, the endplates are reserved and secured between two holding pads to prevent testing artefacts such as friction between these structures. Chemically alterations of discs can be done to promote or recover degenerative degrees. The required eccentric loading compressor will be designed, and customer manufactured. The compressive shearing test will require compression-shear testing system. These tests will be firstly static, and then adapting the static stress damage test results to estimate the fatigue failure load, thus the fatigue loading test can be performed. Digital image correlation may also be introduced to study the morphology and strain distribution of the whole specimen. Digital image correlation technology generally requires those skilled to spray speckle, capture the speckle motion in the loading process by high-speed camera, and calculate the displacement field and strain field by software. The inventor sprays speckle on the entire outer surface of the IVD, and sets up a series of high-speed cameras around the specimen to capture it, so as to obtain and calculate the displacement field and strain field of the whole specimen.

[0059] Clinical assessment

[0060] Finally, clinically assess the application of the model in case-controlled cohort.

[0061] i) Inclusion and exclusion criteria

[0062] For the purposes of this study, 300 patients (aged 18-60 years) with chronic LBP (>3 months history of back pain) have already been recruited with whole spine standing posteroanterior and lateral radiographs and lumbosacral MRI. These patients were recruited as part of a chronic LBP cohort from a specialty adult spine clinic at the Queen Mary Hospital, Hong Kong. These patients all had axial lumbar low back pain with or without symptoms of nerve compression such as radiculopathy, sciatica or lower limb  neurological deficit. All patients with history of previous lumbar spine surgery, tumor, infection or fracture were excluded. Any patient who was unable to undergo radiographs or MRI, such as claustrophobia, resting tremors, pregnant, etc. were also excluded. Written informed consent was obtained from all patients. Ethics was approved by the institutional review board.

[0063] ii) Physiological testing and clinical criteria

[0064] Patients are assessed clinically for location of the back pain and any involvement of one or both lower limbs. Location of the radiculopathy in the lower limbs will be recorded by dermatomal involvement. Lower limb neurological status in terms of numbness pattern and motor deficit in specific myotomes will also be recorded. The severity of pain is measured by the visual analogue scale (VAS) which assesses the patient’s subjective report of the pain from a 0-10 scale. Severity of disability will be assessed by the Oswestry Disability Index (ODI) .

[0065] iii) Imaging assessment criteria

[0066] On radiographs, patients will be assessed for their spino-pelvic alignment parameters. Parameters under study include lumbar lordosis (L1-S1) , pelvic tilt, sacral slope, pelvic incidence, sagittal vertical axis, and any spondylolisthesis and grade of slip on the lateral view; and C7-central sacral vertical line and Cobb angle for scoliosis on the posteroanterior view. On the MRI, the severity of disc degeneration will be assessed by the Pfirmann and Schneidermann scores. Disc herniation will be classified as either disc bulge (base of herniation wider than the apex) , extrusion (base of herniation is narrower than the apex) or sequestration (herniated material detached from main bulk of disc) . The position of the herniation will be assessed on the axial cuts of the MRI and described as anterior, anterolateral, posterior, posterolateral or far lateral.

[0067] Fig. 5 shows a flowchart of the research plan including experiments, computer modelling and clinical validation.

[0068] Step 1: achieving artificial intelligence (AI) assisted segmentation of  degenerative IVD. On the one hand, determining the morphology of degenerative IVD with different level of degeneration, and on the other hand, distinguishing the NP and AF of IVD for the next step of point clouds modeling

[0069] Step 2: redeveloping point clouds model of IVD based on the medical imaging in Step 1. And simulating specific disc biomechanical alterations using point clouds-BSE (point clouds-based computing program) in response to physiologically relevant loading and motion, so as to estimate the severity and progression of degenerative IVD herniation.

[0070] Step 3: conducting biomechanical tests in vitro to examine the association between different morphology of degenerative IVD and altered biomechanical characteristics using material testing system. And Step 2 and 3 verify each other.

[0071] Step 4: clinically evaluating the ability of this system and method in combination with Steps 2 and 3 to estimate the progression of disc degeneration and herniation based on non-invasive radiological scans of patients.

[0072] According to another embodiment of the present invention, a computer-readable storage medium is provided. The computer-readable storage medium comprises a program or executable commands being configured to control a processor to perform the simulation method described above, the implementation principle of which is similar to the above. Instructions or software to control a processor or computer to implement the hardware components and perform the methods as described above are written as computer programs, code segments, instructions or any combination thereof, for individually or collectively instructing or configuring the processor or computer to operate as a machine or special-purpose computer to perform the operations performed by the hardware components and the methods as described above. In one example, the instructions or software include machine code that is directly executed by the processor or computer, such as machine  code produced by a compiler. In another example, the instructions or software include higher-level code that is executed by the processor or computer using an interpreter. Programmers of ordinary skill in the art can readily write the instructions or software based on the block diagrams and the flow charts illustrated in the drawings and the corresponding descriptions in the specification.

[0073] According to another embodiment of the present invention, a system for simulating an intervertebral disc is provided. The system comprises a processor and a storage, wherein, the storage is configured to store executable commands; and the processor is configured to realize the method described above by execute the executable commands.

[0074] The critical invention is a system and method to simulate spine soft tissue (especially the IVD) material properties under loading. This is a system that uses computer simulation to predict how IVD changes due to physiological or pathological stimulations. The system may be used either clinically or non-clinically. In clinical setting, 3D diagnostic images (CT and MRI images) of a patient’s soft tissues, comprising but not limited to IVD, are imported into a computer application that allows surgeons to simulate and compare the biomechanical performance of IVD under physiological or pathological stimulations and to choose different fixation implants, based on the recommendations provided by the simulation results in terms of the optimal types, sizes, and placements of implants to be used. In addition, the mechanical behavior of the diseased disc at all stages can be assessed to further predict the progression of the disease, thus providing targeted preventive treatment methods. In non-clinical settings, the system may be used by R&D engineers to study and improve the designs of orthopedic implants relating to soft tissue. In both clinical and non-clinical applications, simulations are completed using a novel point clouds-based approach to representing soft tissue that was developed by the inventors.

[0075] Embodiments of the present invention have been described in terms of the  preferred embodiment, and it is recognized that equivalents, alternatives, and modifications, aside from those expressly stated, are possible and within the scope of the appending claims.

Claims

1.A method of simulating an intervertebral disc, comprising:automatically segmenting a medical image of a spine based on a pre-trained model so as to obtain different intervertebral disc segmentations; andbiomechanical modelling under multiple motions for the different intervertebral disc segmentations so as to obtain dynamic properties of the intervertebral disc.2.The method of claim 1, wherein the motions comprise a force loading.3.The method of claim 2, wherein, the force loading is simulated with a meshfree method.4.The method of claim 3, wherein, the meshfree method comprises point clouds-based simulation.5.The method of claim 1, wherein, the intervertebral disc segmentations comprise a normal disc, a bulging degenerative disc, a herniated disc and / or a vertebra.6.The method of claim 1, wherein, the intervertebral disc comprises an annulus fibrosis and a nucleus pulposus.7.The method of claim 6, wherein the nucleus pulposus is simulated using a Neo-Hookean model.8.The method of claim 7, wherein a Neo-Hookean strain energy potential  of the Neo-Hookean model is:  where U is the Neo-Hookean strain energy potential; C10 is material constants characterizing the deviatoric deformation of the material; and D1 is material incompressibility parameter; is the first deviatoric strain invariant; Jel is the elastic volume ratio.9.The method of claim 6, wherein the annulus fibrosis is simulated using a Holzapfel-Gasser-Ogden model to represent the nonlinear mechanical behavior of collagen fibers.10.The method of claim 9, wherein a Holzapfel-Gasser-Ogden strain energy potential of the Holzapfel-Gasser-Ogden model is as follows: wherein, andU is the Holzapfel-Gasser-Ogden strain energy potential; C10 is material constants characterizing the deviatoric deformation of the material; D1 is material incompressibility parameter; k1 is a dimensionless parameter; k2 is a stress-like parameter to be determined from mechanical tests of the tissue; κ is structure parameters to be determined from histological data of the tissue; N is the number of families of fibers; is the first deviatoric strain invariant; Jel is the elastic volume ratio.11.The method of claim 10, wherein, when biomechanical modelling under multiple motions for the different intervertebral disc segmentations, considering a fiber mechanics model expressed as follows: wherein, σf and εf correspondingly represent fiber stress and strain, representing elasticity and creep, df is a damage variable, 0≤df≤1, 0 representing no damage to the fiber, and 1 representing completely tear.12.The method of claim 1, wherein the pre-trained model is a deep learning model, which are trained by using the samples of medical images of spines.13.The method of claim 12, wherein the model is trained by the following steps:Collecting the medical images of spines;Segmenting, masking and labeling intervertebral disc, nucleus, and spinal canal of the images, to obtain the samples of medical images of spines;training the deep learning model by using the samples of medical images of spines.14.A computer-readable storage medium comprising a program, the program is executed to perform the method of any one of claims 1-13.15.A system for simulating an intervertebral disc, comprising a processor and a storage, wherein, the processor is configured to perform the method of any one of claims 1-13 by execute the executable commands in the storage.

Images