Spatiotemporal Manifold Prediction Model for Anterior Vertebral Body Growth Modulation Surgery in Idiopathic Scoliosis

Spatiotemporal Manifold Prediction Model for Anterior Vertebral Body Growth Modulation Surgery in Idiopathic Scoliosis

William Mandel 1MedICAL, Polytechnique Montreal, Montreal, QC, Canada 1    Olivier Turcot 2Sainte-Justine Hospital Research Center, Montreal, QC, Canada 2    Dejan Knez 3University of Ljubljana, Faculty of Electrical Engineering, Slovenia 13   
Stefan Parent
2Sainte-Justine Hospital Research Center, Montreal, QC, Canada 2
   Samuel Kadoury 1MedICAL, Polytechnique Montreal, Montreal, QC, Canada 12Sainte-Justine Hospital Research Center, Montreal, QC, Canada 2
thanks: Supported by the Canada Research Chairs and NSERC Discovery Grants.

Anterior Vertebral Body Growth Modulation (AVBGM) is a minimally invasive surgical technique that gradually corrects spine deformities while preserving lumbar motion. However the selection of potential surgical patients is currently based on clinical judgment and would be facilitated by the identification of patients responding to AVBGM prior to surgery. We introduce a statistical framework for predicting the surgical outcomes following AVBGM in adolescents with idiopathic scoliosis. A discriminant manifold is first constructed to maximize the separation between responsive and non-responsive groups of patients treated with AVBGM for scoliosis. The model then uses subject-specific correction trajectories based on articulated transformations in order to map spine correction profiles to a group-average piecewise-geodesic path. Spine correction trajectories are described in a piecewise-geodesic fashion to account for varying times at follow-up exams, regressing the curve via a quadratic optimization process. To predict the evolution of correction, a baseline reconstruction is projected onto the manifold, from which a spatiotemporal regression model is built from parallel transport curves inferred from neighboring exemplars. The model was trained on 438 reconstructions and tested on 56 subjects using 3D spine reconstructions from follow-up exams, with the probabilistic framework yielding accurate results with differences of in main curve angulation, and generating models similar to biomechanical simulations.

1 Introduction

Spinal morphology and more particularly 3D morphometric parameters, have demonstrated significant potential in assessing the risk of spinal disease progression. For spinal deformities such as adolescent idiopathic scoliosis (AIS), personalized 3D reconstructions generated from radiographs allows surgeons to assess the severity and decide on efficient treatment options. A recently introduced minimally invasive surgical technique called Anterior Vertebral Body Growth Modulation (AVBGM) consists of instrumenting the spine with traditional vertebral implants to link a segment of vertebrae together by a flexible polypropylene cable applied to the spine anteriorly. The fusion-less technique applies compressive forces on the convex side of the spinal curve, thereby modulating the distribution of pressure on the vertebral growth plates. In combination with natural bone growth, this allows to retain spine flexibility [skaggs2014classification]. While this new surgical technique showed promising results for skeletally immature patients [crawford2010growth], difficulties were reported to predict short and long-term post-operative correction [samdani2015anterior]. Biomechanical models were shown to reproduce surgical outcomes, but are not adapted for real-time surgical applications [cobetto20183d].

A recent study evaluated differences in hand-crafted 3D parameters in progressive AIS groups using images from the patient’s first visit [nault2014three] to predict progression, by manually selecting the best features that can characterize the intrinsic nature of 3D spines. On the other hand, dimensionally-reduced growth trajectories of various anatomical sites have been investigated in neurodevelopment studies for newborns, based on geodesic shape regression to compute the diffeomorphisms based on image time series of a population [singh2013hierarchical]. These regression models were also used to estimate spatio-temporal evolution of the cerebral cortex, by automatically identifying the points of interest and inertia between the first and follow-up images based on non-rigid transformations [fishbaugh2014geodesic]. The concept of parallel transport curves in the tangent space from low-dimensional manifolds proposed by Schiratti et al. [schiratti2015learning] was used to analyze shape morphology [kadoury20173] and adapted for radiotherapy response [chevallier2017learning], but lacks the capability to predict correction from applied forces following surgery.

This paper presents a prediction model for patient response to AVBGM from pre-operative 3D spine models reconstructed from biplanar X-ray images (Fig. 1). The method first trains a piecewise-geodesic manifold using a collection of pre-operative and longitudinal 3D reconstructions of the spine acquired during follow-up evaluations of patients treated with AVBGM for AIS. A discriminant adjacency matrix is constructed to separate responding and non-responding patients. During testing, an unseen baseline spine model is projected onto the manifold, where a piecewise-geodesic curve describing spatiotemporal evolution is regressed using discrete approximations, from which the curvature evolution is inferred, yielding a prediction of the intervertebral displacements and shape morphology describing deformation correction. The main contribution of this paper is the introduction of a piecewise-geodesic transport curve in the tangent space from low-dimensional samples designed for the correction of spinal deformities, where a new time-warping function controlling the rate of correction is obtained from clinical parameters.

2 Method

2.1 Discriminant Embedding of Longitudinal Spine Models

Figure 1: Proposed prediction framework for spine surgery outcomes. In the training phase, a dataset of spine models are embedded in a spatio-temporal manifold , into responsive (R) or non-responsive (NR) groups. During testing, an unseen baseline 3D spine reconstruction is projected on using based on Nadaraya-Watson kernels. The closest samples to the projected point x are selected to regress the spatiotemporal curve used for predicting the correction due with AVBGM.

A sample spine reconstruction is represented by , modelling a series of vertebral shapes. For each vertebra, template-based models are obtained where vertex coordinates have one-to-one correspondences between samples. In addition to the mesh-based representation, each model possesses a list of annotated landmarks to compute local intervertebral rigid transformations, such that , with a rigid inter-vertebral transform. Hence, the overall shape of the spine is described as a vector of sequential registrations assigned to each vertebral level, whereby considering the ensemble of transformations, we obtain a combination of previous transforms:


using recursive compositions. The feature array dictates the location and rotation of the object constellation, while procuring the morphology of the vertebrae S. The model is deformed by applying displacements to the inter-vertebral parameters. By extending this to the entire absolute vector representing the spine model, this then achieves a global deformation. In this case, registrations are described in the reference coordinate system of the lower vertebra, corresponding to it’s principal axes of the cuneiform shape with the origin positioned at the center of mass of the vertebra. The rigid transformations are the combination of a rotation matrix and a translation vector . We formulate the rigid transformation of a vertebra mesh as where . Composition is given by .

We propose to embed a collection of non-responsive (NR) and (2) responsive (R) patients to AVBGM which will offer a maximal separation between the classes, by using a discriminant graph-embedding. Here, labelled points defined in are embedded in the low-dimensional manifold , where describes the label (NR or R) and defines the time of follow-up. We assume that for the sampled data, an underlying manifold of the high-dimensional data exists such that defined in . We rely on the assumption that a locally linear mapping exists, where local neighbourhoods are defined as tangent planes estimated with and , describing the paired distances between linked neighbours . Hence, the relationship can be established as .

Because the discriminant manifold structure in requires to maintain the local structure of the underlying data, a undirected similarity graph is built, where each node V are connected to each other with edges that are weighted with the graph W. The overall structure of is therefore defined with for feature vectors belonging to the same class and , which separate features from both classes. During the embedding of the discriminant locally linear latent manifold, data samples are divided between and .

2.2 Piecewise-Geodesic Spatiotemporal Manifold

Once sample points are in manifold space, the objective is to regress a regular and smooth piecewise-geodesic curve that accurately fits the embedded data describing the spatiotemporal correction following AVBGM within a 2 year period. For each sample data , the closest individuals demonstrating similar baseline features are identified from the embedded data, creating neighborhoods with measurements at different time points, thus creating a low-dimensional Riemannian manifold where data points , with denoting a particular individual, the time-point measurement and the pre-operative model. By assuming the manifold domain is complete and piecewise-geodesic curves are defined for each time trajectories, time-labelled data can be regressed continuously in , thereby creating smooth curves in time intervals described by samples in .

However, due to the fact the representation of the continuous curve is a variational problem of infinite dimensional space, the implementation follows a discretization process which is derived from the procedure in [boumal2011discrete], such that:


This minimization process simplifies the problem to a quadratic optimization, solved with LU decomposition. The piecewise nature is represented by the term , defined as samples along . The 1 component of Eq.(2.2) is a penalty term to minimize the geodesic distance between samples and the regressed curve, where are weight variables based on sample distances. This helps regress a curve that will lie close to , shifted by in order to have the initial reconstructions co-registered. The 2 term represents the velocity of the curve (defined by , approximating ), minimizing the distance of the 1 derivative of . By minimizing the value of the curve’s first derivatives, this prohibits any discontinuities or rapid transitions of the curve’s direction, and is modulated by . Finally, an acceleration penalty term (defined by ) focuses on the 2 derivative of with respect to by minimizing the norm. The acceleration is modulated by . Estimates for and (weighted by , respectively), are generated using geometric finite differences. These estimates dictates the forward and backward step-size on the regressed curve, leading to directional vectors in as shown in [boumal2011discrete]. In order to minimize , a non-linear conjugate gradient technique defined in the low-dimensional space is used, thus avoiding convergence and speed issues. The regressed curve is therefore defined for all time points, originating at . The curve creates a group average of spatiotemporal transformations based on individual correction trajectories.

2.3 Prediction of Spine Correction

Finally, to predict the evolution of spine correction from an unseen pre-operative spine model, we use the geodesic curve modelling the spatiotemporal changes of the spine, where each point is associated to a speed vector v defined with a tangent plane on the manifold such that .

Based on Riemannian theory, an exponential mapping function at x with velocity v can be defined from the geodesics such that . Using this concept, parallel transport curves defined in can help define a series of time-index vectors along as proposed by [schiratti2015learning]. The collection of parallel transport curves allows to generate an average trajectory in ambient space , describing the spine changes due to the corrective forces of tethering. The general goal is to begin the process at the pre-operative sample, and navigate the piecewise-geodesic curve describing correction evolution in time, where one can extract the appearance at any point (time) in using the exponential mapping. For implementation purposes, the parallel transport curve are constrained within a smooth tubular boundary perpendicular to the curve (from an ICA) to generate the spatiotemporal evolution in the coordinate system of the pre-operative model.

Hence, given the manifold at time with v defined in the tangent plane and the regressed piecewise-geodesic curve , the parallel curve is obtained as:


Therefore by repeating this mapping for manifold points seen as samples of individual progression trajectories along , an evolution model can be generated. Whenever a new sample is embedded, new samples points along , denoted as can be generated parallel to the regressed piecewise curve in , capturing the spatiotemporal changes in correction.

A time warp function allowing to vary along the geodesic curve is described as . Here, we propose to incorporate a personalized acceleration factor based on the spine maturity and flexibility derived from the spine bending radiographs and Risser grade. A coefficient describing the change in Cobb angle between poses, and modulated by the Risser grade . This coefficient regulates the rate of correction based on the neighbouring samples. Finally, to take under account the relative differences between the group-wise samples and the query model once mapped onto the regressed curve, a time-shift parameter is incorporated in the warp function.

For spine correction evolution, displacement vectors are obtained by a PCA of the hyperplane crossing in manifold [schiratti2015learning]. Hence, for any query sample which represents the mapped pre-operative 3D reconstruction (prior to surgery), the predicted model at time can be regressed from the piecewise-geodesic curve generated from embedded samples x in such that:


which yields a predicted post-operative model in high-dimensional space , and a zero-mean Gaussian distribution. The generated model offers a complete constellation of inter-connected vertebral models composing the spine shape S, at first-erect (FE), 1 or 2-year visits, including landmarks on vertebral endplates and pedicle extremities, which can be used to capture the local shape morphology with the correction process.

3 Experiments

The discriminant manifold was trained from a database of 3D spine reconstructions generated from biplanar images [Humbert09], originating from patients demonstrating several types of deformities with immediate follow-up (FE), 1 year and 2 year visits. Patients were recruited from a single center prospective study, with the inclusion criteria being evaluated by an orthopaedic surgeon and a main curvature angle between 30 and 60. Patients were divided in two groups, with the first group composed of 94 responsive patients showing a reduction in Cobb angle over or equal to 10 between the FE and follow-up visit. The second group was composed of 37 non-responsive (NP) patients with a reduction of less than 10. Each vertebra model of the spine were annotated with 4 pedicle tips and 2 center points placed on the vertebral endplates, and validated by an experienced radiologist. These expert-selected landmarks were used to establish the local coordinate system for each vertebra, describing the orientation and location (known pose), and used as control points to warp triangulated shape models generated from CT images of a cadaveric spine.

We evaluated the geometrical accuracy of the predictive manifold for 56 unseen surgical patients with AVBGM (mean age , average main Cobb angle on the frontal plane at the first visit was ), with predictions at (FE), and months. For the predicted models, we evaluated the 3D root-mean-square difference of the vertebral landmarks generated, the Dice coefficients of the vertebral shapes and in the main Cobb angle. The results are shown in Table 1. Results were confronted to other techniques such as biomechanical simulations performed on each subject using finite element modelling with ex-vivo parameters [cobetto20183d], a locally linear latent variable model [park2015bayesian] and a deep auto-encoder network [thong2016three]. Fig. 2 shows a sample prediction result for an 11 y.o. patients at FE, 12 and 24-months for a patient with right-thoracic deformity, which are more common in the scoliotic population. Results from the predicted geometrical models show the regressed spatio-temporal geodesic curve yields anatomically coherent structures, with accurate local vertebral morphology.

To evaluate robustness with respect to varying instrumented levels, we measured the accuracy of the predicted models for tethering between 4 and 8 vertebrae at 2 yrs, ranging from thoracic to lumbar regions. Fig. 2(b) shows the improvement of the spatiotemporal geodesic curve in comparison to traditional biomechanical models, particularly when the number of levels are higher.

4 Conclusion

In this paper, we proposed an accurate predictive model of spine morphology and Cobb angle correction obtained at the first-erect, 1-year and 2-year visits, following anterior vertebral body growth modulation. The piecewise-geodesic curve capturing spatio-temporal changes could be used for patient selection of AVBGM as a decision-sharing tool prior to surgery. Our approach is based on smooth and regular trajectories embedded in a discriminant manifold, which enable an efficient navigation on a low-dimensional domain trained from operative cases, yielding results similar to actual surgical outcomes. Future work will include a multi-center evaluation before it can be used in clinical practice.

FE visit 1-year visit 2-year visit
3D RMS Dice Cobb 3D RMS Dice Cobb 3D RMS Dice Cobb
Biomec. sim 3.31.1 853.4 2.80.8 3.61.2 843.6 3.20.9 4.12.3 823.9 3.61.0
LL-LVM [park2015bayesian] 3.61.4 834.0 3.81.5 4.73.3 794.4 5.52.6 6.64.4 715.9 7.03.9
Deep AE [thong2016three] 4.11.5 804.4 5.12.7 5.01.9 774.9 5.83.0 6.34.6 725.7 6.64.2
Proposed 2.40.8 922.7 1.80.5 2.90.9 902.8 2.00.7 3.21.3 873.1 2.10.6
Table 1: 3D RMS errors (mm), Dice (%) and Cobb angles () for the proposed method, and compared with biomechanical simulations, locally linear latent variable models (LL-LVM) and deep auto-encoders (AE). Predictions are evaluated at FE, 1 and 2-yrs.
(a) (b)
Figure 2: (a) Comparison in actual and predicted Cobb angles in a 11 y.o. patient at the first-erect visit, at 1-yr and at 2-yrs postop. Top row depicts the actual X-rays, while the bottom row presents the predicted 3D spine geometry. (b) Errors with 5 different tethering levels, comparing results with biomechanical simulations at 2 yrs.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description