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

Abstract

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 $2.1\pm 0.6^{\circ }$ in main curve angulation, and generating models similar to biomechanical simulations.

Supported by the Canada Research Chairs and NSERC Discovery Grants.

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 [1]. While this new surgical technique showed promising results for skeletally immature patients [2], difficulties were reported to predict short and long-term post-operative correction [3]. Biomechanical models were shown to reproduce surgical outcomes, but are not adapted for real-time surgical applications [4].

A recent study evaluated differences in hand-crafted 3D parameters in progressive AIS groups using images from the patient’s first visit [5] 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 [6]. 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 [7]. The concept of parallel transport curves in the tangent space from low-dimensional manifolds proposed by Schiratti et al. [8] was used to analyze shape morphology [9] and adapted for radiotherapy response [10], 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

A sample spine reconstruction is represented by $\mathbf S = \{\mathbf {s}_{1},\ldots ,\mathbf {s}_{m}\}$, modelling a series of $m=17$ vertebral shapes. For each $\mathbf {s}_{i}$ 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 $\mathbf {s}_{i}$ possesses a list of annotated landmarks to compute local intervertebral rigid transformations, such that $A = [T_{1}, T_{2}, \ldots , T_{m}]$, with $T_i=\{R,t\}$ 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:

$$\begin{aligned} \mathbf y _{\text {i}} = [T_{1}, \mathbf {s}_{1}; T_{1} \circ T_{2}, \mathbf {s}_{2},\ldots , T_{1} \circ T_{2} \circ \ldots \circ T_{m},\mathbf {s}_{m}] \end{aligned}$$

(1)

using recursive compositions. The feature array $\mathbf y _{\text {i}}$ dictates the location and rotation of the object constellation, while procuring the morphology of the vertebrae $\mathbf 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 R and a translation vector t. We formulate the rigid transformation $T = \{R, t\}$ of a vertebra mesh $\mathbf {s}_{i}$ as $y = Rx+t$ where $x, y, t \in \mathfrak {R}^{3}$. Composition is given by $T_{1}\circ T_{2}=\{R_{1}R_{2},R_{1}t_{2}+t_{1}\}$.

Fig. 1.

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

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, n labelled points $\mathbb {Y}=\{(\mathbf y _i,l_i,t_i)\}_{i=1}^{n}$ defined in $\mathbb {R}^{D}$ are embedded in the low-dimensional manifold $\mathcal {M}$, where $l_i$ describes the label (NR or R) and $t_i$ defines the time of follow-up. We assume that for the sampled data, an underlying manifold of the high-dimensional data exists such that $\mathbb {X}=\{(\mathbf x _i,l_i,t_i)\}_{i=1}^{n}$ defined in $\mathbb {R}^{d}$. We rely on the assumption that a locally linear mapping $\mathbf M _i \in \mathbb {R}^{D \times d}$ exists, where local neighbourhoods are defined as tangent planes estimated with ${ \mathbf y _j - \mathbf y _i}$ and ${ \mathbf x _j - \mathbf x _i}$, describing the paired distances between linked neighbours i, j. Hence, the relationship can be established as $\mathbf y _j - \mathbf y _i \approx \mathbf M _i ( \mathbf x _j - \mathbf x _i )$.

Because the discriminant manifold structure in $\mathbb {R}^{d}$ requires to maintain the local structure of the underlying data, a undirected similarity graph $\mathcal {G}=(\varvec{{V}},\varvec{{W}})$ is built, where each node $\varvec{{V}}$ are connected to each other with edges that are weighted with the graph $\varvec{{W}}$. The overall structure of $\mathcal {M}$ is therefore defined with ${\varvec{W}}_w$ for feature vectors belonging to the same class and ${\varvec{W}}_b$, which separate features from both classes. During the embedding of the discriminant locally linear latent manifold, data samples are divided between ${\varvec{W}}_w$ and ${\varvec{W}}_b$.

2.2 Piecewise-Geodesic Spatiotemporal Manifold

Once sample points $\mathbf x _i$ are in manifold space, the objective is to regress a regular and smooth piecewise-geodesic curve $\gamma : [t_1,t_N]$ that accurately fits the embedded data describing the spatiotemporal correction following AVBGM within a 2 year period. For each sample data $\mathbf x _i$, the K closest individuals demonstrating similar baseline features are identified from the embedded data, creating neighborhoods $\mathcal {N}(\mathbf x _q)$ with measurements at different time points, thus creating a low-dimensional Riemannian manifold where data points $\mathbf x _{i,j}$, with i denoting a particular individual, j the time-point measurement and $j=0$ 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 $\mathbb {R}^{D}$, thereby creating smooth curves in time intervals described by samples in $\mathbb {R}^{d}$.

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 [11], such that:

$$\begin{aligned} E(\gamma ) = \dfrac{1}{K_d} \sum _{i=1}^{K_d}&\sum _{j=0}^{t_N} w_i \Vert \gamma (t_{i,j}) - (\mathbf x _{i,j} - (\mathbf x _{i,0} - \mathbf x _q)) \Vert ^2 \nonumber \\&+ \dfrac{\lambda }{2} \sum _{i=1}^{K_d} \alpha _i \Vert v_i\Vert ^2 + \dfrac{\mu }{2} \sum _{i=1}^{K_d} \beta _i \Vert a_i \Vert ^2. \end{aligned}$$

(2)

This minimization process simplifies the problem to a quadratic optimization, solved with LU decomposition. The piecewise nature is represented by the term $K_d \in \mathcal {N}(\mathbf x _q)$, defined as samples along $\gamma$. The 1$^{st}$ component of Eq. (2) is a penalty term to minimize the geodesic distance between samples $\mathbf x _{i,j}$ and the regressed curve, where $w_i$ are weight variables based on sample distances. This helps regress a curve that will lie close to $\mathbf x _{i,j}$, shifted by $\mathbf x _q$ in order to have the initial reconstructions co-registered. The 2$^{nd}$ term represents the velocity of the curve (defined by $v_i$, approximating $\dot{\gamma }(t_i)$), minimizing the $L_2$ distance of the 1$^{st}$ derivative of $\gamma$. 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 $\alpha _i$. Finally, an acceleration penalty term (defined by $a_i$) focuses on the 2$^{nd}$ derivative of $\gamma$ with respect to $t_i$ by minimizing the $L_2$ norm. The acceleration is modulated by $\beta _i$. Estimates for $v_i$ and $a_i$ (weighted by $\{ \lambda ,\mu \}$, 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 $\mathcal {M}$ as shown in [11]. In order to minimize $E(\gamma )$, a non-linear conjugate gradient technique defined in the low-dimensional space $\mathbb {R}^{d}$ is used, thus avoiding convergence and speed issues. The regressed curve $\gamma$ is therefore defined for all time points, originating at $t_0$. 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 $\gamma : \mathbb {R}^{D} \rightarrow \mathcal {M}$ modelling the spatiotemporal changes of the spine, where each point $\mathbf x \in \mathcal {M}$ is associated to a speed vector $\mathbf v$ defined with a tangent plane on the manifold such that $\mathbf v \in \text {T}_\mathbf{x } \mathcal {M}$.

Based on Riemannian theory, an exponential mapping function at $\mathbf x$ with velocity $\mathbf v$ can be defined from the geodesics such that $\text {e}_\mathbf{x }^{\mathcal {M} }(\mathbf v )$. Using this concept, parallel transport curves defined in $\text {T}_\mathbf{x }$ can help define a series of time-index vectors along $\gamma$ as proposed by [8]. The collection of parallel transport curves allows to generate an average trajectory in ambient space $\mathbb {R}^{D}$, 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 $\mathbb {R}^{D}$ 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 $t_0$ with $\mathbf v$ defined in the tangent plane and the regressed piecewise-geodesic curve $\gamma$, the parallel curve is obtained as:

$$\begin{aligned} \eta ^\mathbf{v } (\gamma ,s)= \text {e}_{ \gamma (s)}^{\mathcal {M}} (\mathbf x _{\gamma ,t_0, s}(\mathbf v )), \,\,\, s \in \mathbb {R}^{d}. \end{aligned}$$

(3)

Therefore by repeating this mapping for manifold points seen as samples of individual progression trajectories along $\gamma (s)$, an evolution model can be generated. Whenever a new sample is embedded, new samples points along $\gamma (s)$, denoted as $\eta ^\mathbf{v } (\gamma ,\cdot )$ can be generated parallel to the regressed piecewise curve in $\mathcal {M}$, capturing the spatiotemporal changes in correction.

A time warp function allowing s to vary along the geodesic curve is described as $\phi _i(t)= \theta _i (t-t_0-\tau _i) +t_0$. 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 $\theta _i = C_i \times R_i$ describing the change in Cobb angle $C_i$ between poses, and modulated by the Risser grade $R_i$. This coefficient regulates the rate of correction based on the K 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 $\tau _i$ is incorporated in the warp function.

For spine correction evolution, displacement vectors $\mathbf v _i$ are obtained by a PCA of the hyperplane crossing $\text {T}_\mathbf{x _i} \mathcal {M}$ in manifold $\mathcal {M}$ [8]. Hence, for any query sample $\mathbf x _q$ which represents the mapped pre-operative 3D reconstruction (prior to surgery), the predicted model at time $t_k$ can be regressed from the piecewise-geodesic curve generated from embedded samples $\mathbf x$ in $\mathcal {N}(\mathbf x _q)$ such that:

$$\begin{aligned} \mathbf y _{q,t_k} = \eta ^\mathbf{v _q} (\gamma ,\phi _i(t_{k})) + \epsilon _{q,t_k} \end{aligned}$$

(4)

which yields a predicted post-operative model $\mathbf y _{q,t_k}$ in high-dimensional space $\mathbb {R}^{D}$, and $\epsilon _{q,t_k}$ a zero-mean Gaussian distribution. The generated model offers a complete constellation of inter-connected vertebral models composing the spine shape $\mathbf 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 438 3D spine reconstructions generated from biplanar images [12], originating from 131 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$^{\circ }$ and 60$^{\circ }$. 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$^{\circ }$ 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$^{\circ }$. 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 $12 \pm 3$, average main Cobb angle on the frontal plane at the first visit was $47^{\circ } \pm 10^{\circ }$), with predictions at $t=0$ (FE), $t=12$ and $t=24$ 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 [4], a locally linear latent variable model [13] and a deep auto-encoder network [14]. Figure 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. Figure 2(b) shows the improvement of the spatiotemporal geodesic curve in comparison to traditional biomechanical models, particularly when the number of levels are higher.

Table 1. 3D RMS errors (mm), Dice (%) and Cobb angles ($^{\circ }$) 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.

Fig. 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.

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.