Patient Specific Biomechanics Are Clinically Significant In Accurate Computer Aided Surgical Image Guidance

Our patient specific AR IG method uses a mesh free physical simulation. Mesh free approaches are suited to AR IG because they do not require re-meshing after topology changes (such as cutting) that occur during a resection and have been previously implemented in SOFA [17]. Our MRE simulation uses the following canonical semi-discrete formulation:

where $\mathbf{q}$ is the nodal displacement vector, $\mathbf{\dot{q}}$ is the nodal velocity vector, $\mathbf{\ddot{q}}$ is the nodal acceleration vector, $\mathbf{M}$ is the mass matrix, $\mathbf{f_{int}}$ is the internal force vector, and $\mathbf{f_{ext}}$ is the external force vector.

We choose a computationally efficient discretization of this system by adopting the Implicit Euler method [18] which will create a system of linear equations. This is because our cohort MRE data encodes linear mechanics and because linear approximations are used to achieve video rate performance in AR IG. The system of equations is fed into the Conjugate Gradient (CG) descent method that will ultimately solve for $\delta{\mathbf{\dot{q}}}$ for each time step. The implicit Euler method is adopted as follows:

where $h$ is the time step (this is in milliseconds for our study), $\mathbf{K}$ is the Jacobian of first order derivatives of $\mathbf{f_{int}}$ with respect to $\mathbf{\dot{q}}$ (also known as the stiffness matrix), and $\mathbf{C}$ is the Jacobian of first order derivatives of $\mathbf{f_{int}}$ with respect to $\mathbf{q}$ (also known as the damping matrix). Note that $K$, $C$ and $M$ may be precomputed offline and that this corresponds to “building the simulator model” from processed preoperative MRI and MRE in fig. 3.

The next step of fig. 3, Physical Simulation, is online. Equation 2 is solved numerically via the following mappings:

Such that $A\mathbf{x}=\mathbf{b}$, where $\mathbf{x}$ contains the vector encoding of surgical landmark positions to resolve for the current time step. Once $\delta\mathbf{\dot{q}}$ is solved for each node in the system, calculating the change in global position, $\delta\mathbf{q}$, is simply done by integration.

Our MRE method can improve AR IG accuracy with no computational overhead since it can be shown that online computation method of $\delta\mathbf{\dot{q}}$ is the bottleneck. First, there is an initial guess: $\mathbf{x_{0}}$

(When $n=0$, $\mathbf{p_{0}}$ is initially set to $\mathbf{r_{0}}$)

The steps above illustrate a portion of one iteration of the conjugate gradient descent method, with each subsequent $\mathbf{x_{n}}$ representing a closer approximation of $\mathbf{x}$. The linear solver continues to iterate until the loop is terminated upon reaching preset error bounds (in the example given, this happens when $n=N$ which is 200 in our study). Since MRE data is used to construct $A$ it does not exacerbate the CG bottleneck.

We create a SOFA [19] Physics Model “node” that can be integrated into any SOFA AR IG system. The node is composed of a stiffness matrix ($\mathbf{K}$) created from MRE data, geometry created from MRI data, a mass matrix ($\mathbf{M}$) computed using MRE and MRI data and a displacement vector $\mathbf{q}$ that is initialized from MRI data and updated by the physical simulation during AR IG.

Patient liver models are built using the procedure depicted in fig. 3 (gray region). Firstly liver data from our cohort MRI and MRE scans is segmented using the tool described in fig. 2. Next surface geometry is created from MRI data using a modified Delaunay triangulation [20]. Liver biomechanics are then extracted from MRE data. We first convert voxel intensity to elastic modulus, then sub-sample this for real time AR IG using a 3D Voronoi decomposition of the volume to choose DOFs. After $\mathbf{M}$ is computed, $\mathbf{K}$,$\mathbf{M}$, and $\mathbf{q}$ are assembled into the SOFA node.

Fig. 3 shows the major steps of a retraction experiment using our MRE AR IG method. We compare blood vessel positions from our method with a traditional AR IG method. Our focus is on clinical significance of our landmark placement method and our experimental pipeline differs from the canonical AR IG pipeline in fig. 1. We emulate AR steps prior to physical simulation (fig. 3 green region) because accurate real time implementations of all major AR steps remain open research problems. Laparoscopic video input is noisy and many accuracy trade offs are made in the 3D reconstruction and registration steps that would make analysis of our Physical simulation method difficult. Liver retraction (hoisting) video input is therefore emulated using a virtual implementation of the 10mm diameter Nathanson retractor shown hoisting the liver in fig. 3.

The traditional AR IG approach is an implementation of Plantafeve’s state of the art model [6] that uses a traditional “atlas” reference stiffness from Lee’s study [16]. Traditional AR IG can be directly compared with our data driven simulations by using the same emulated retraction input as our MRE method.

Liver retractions were simulated using the same cohort described in section II-A. We assumed the hoisting force was equal to the gravitational force of the liver and was coplanar on the transverse plane, with the Nathanson retractor force applied to the lower inferior right lobe. We further assume a linear elastic resistance from the abdomen and that friction is negligible. We rejected 17 of our cohort scans where it was not possible to simulate a retraction because of bad geometry, leaving 103 cases of vessel placement compared.

Standard AR IG validation approach is to measure simulation registration error with a phantom (tissue proxy) or a theoretical model [21, 22]. Because MRE validation has been reported using theoretical models [23] we choose a theoretical validation. We simulate a cantilever beam under distributed load and check for convergence of our simulation method with an established theoretical model of the system and our mesh free model.

In addition to testing our simulation method, we use a conventional Finite Element Analysis (FEA) simulation method to model beam deflection as a convergence baseline. This helps better understand the accuracy and speed trade off of our real time mesh free model.

Major steps of our validation study are shown in fig. 4. Our simulation method and the traditional method are implemented using the SOFA framework. The theoretical model is a closed form Euler-Bernoulli model and beam deflection can be computed directly. All three beams model a silicone phantom since the phantom stiffness can be reasonably linearized for a validation model [24], our MRE scans have linear elastic modulus and beam models have been used for SOFA mesh free validation in the literature [17]. We choose a high resolution FEA method with 1.64mm³ element geometry as the baseline.

We use a shape function to compare our mesh free beam to the theoretical beam as they have different coordinate systems. The mapping result can be seen in the “Phantom Physics Model (Our Method)” of fig. 4. The green line is the Voronoi to Euler-Bernoulli coordinate mapping. We then use registration distance to compute the convergence error as is common in the literature [21, 22]. Details of the shape function can be found in [17]. The procedure for the FEA baseline is similar except FEA DOFs coincident with the green line of the “Phantom Physics Model (Traditional Method)” in fig. 4 can be directly compared with the theoretical model. The Euler-Bernoulli beam equation reads:

where $w(x)$ is displacement of the central beam axis in the load direction, $E$ is Young’s modulus, $L$ is the long edge and $I$ is the second moment of beam cross sectional area.

For validation simulation geometry, a hyperrectangle beam with idealized axisymmetric AAA geometry was discretized with dimensions l=50mm, w=10mm, h=10mm. $E$ was the 12kPa Young’s modulus of a well characterized 10% silicone phantom from Clear Ballistics LLC. We set loading and convergence criteria following Plantefeve’s patient specific modeling approach for liver vessel tracking [6].