[ad_1]

Twelve right-handed medical students (four female and eight male) and two professional surgeons (both male) were recruited for this study. We labelled medical students with letters A to L, and surgeons with”SA” and “SB” labels (referring to surgeon A and B, respectively). Only three subjects (A, C and D) repeated the trials (two months after the first trial). Subjects that repeated the trials have a numeral in the label indicating the trial order (e.g. “A2” means the second trial of subject A). None of the student participants had any prior experience in surgical cutting tasks. The study was approved by the University of Edinburgh, School of Informatics, Informatics Ethics panel. All participants provided written informed consent to participate in this study.

The participants were asked to perform a series of 6 elliptical excisions on the phantom using the sensorised cutting tool (Fig. 1d). Before each trial, a new blade (Swann-Morton No. 10) was mounted to the cutting tool. After receiving the task instructions, participants were familiarized with the experimental setup, cutting tool ergonomics, phantom mechanical properties, etc. Next, each subject was asked to rehearse the described task using a dedicated sacrificial phantom. During the trials, the cutting forces that act on the blade in the direction of cutting were recorded at a fixed frequency of 30 Hz. Finally, at the end of the trials, each participant was asked to complete a post-study questionnaire.

Data measurement

Each participant performed six elliptical excisions as a part of the task, yielding 12 force profiles per trial (each excision consists of upper and lower cuts). The recorded profiles were time-aligned and cropped to a fixed duration of 120 samples or 4 seconds (at a sampling rate of 30 Hz). Finally, the samples were normalized to the maximum force value in the entire dataset.

Given the normalized force profiles f(t) (Supplementary Fig. 4a), the virtual displacement profiles x(t) (Supplementary Fig. 4b) were obtained by solving the differential equation for the Maxwell model, equation (1), as follows:

$$x(t)=\frac{f(t)}{E}+\frac{1}{\eta }\int\nolimits_{0}^{T}f(t)\,{{{{{{{\rm{d}}}}}}}}t$$

(1)

where f(t) is the corresponding force profile, T is the duration of the force profile, η = 0.5 N s cm−1 and E = 1 N cm−1 are Maxwell model’s damping and spring coefficients, respectively.

The corresponding virtual velocity profiles \(\dot{x}(t)\) (Supplementary Fig. 4c) were obtained by approximating the time derivative of x(t) using the finite difference method with a step size dt = 0.033.

Elliptical excision force model

The collected incision force profiles show temporal features that can characterize the cutting behaviour. For example, the characteristic dip in an incision force profile (Fig. 1b) might reflect a dynamic change in the configuration of the blade, tissue tensioning applied by the non-dominant hand, or both. Here, we propose a generative model that captures these subject-specific temporal features in the force profiles and enables the disentanglement of skill from incision force analysis.

Figure 1 a shows the approximate model of the task of cutting a viscoelastic phantom as a continuous blade’s movement through a Maxwell body. In the context of this approximation, the Maxwell model29 relates the actual incision force f(t) to a “virtual” velocity of the blade \(\dot{x}(t)\), as follows:

$$\eta \dot{x}(t)=f(t)+\frac{\eta }{E}\dot{f}(t)$$

(2)

where \(\dot{f}(t)\) is the time derivative of the force, and η and E are the Maxwell model’s damping and spring coefficients, respectively.

By taking the Laplace transform of equation (2) and rearranging the terms, we obtain the transfer function G(s), which relates a virtual blade’s displacement X(s) and the actual force F(s), as follows:

$$G(s)=\frac{F(s)}{X(s)}=\frac{\eta s}{\frac{\eta }{E}s+1}$$

(3)

The above transfer function indicates that the model exhibits high-pass characteristics in the force response to the displacement input. This predicts an exponential decay of force with a time constant \(\frac{\eta }{E}\), as a response to a unit step displacement. Importantly, this also predicts a step-like response in the force to a ramp-like displacement input, and therefore, the observed cutting force profiles can be described as a response to a continuous virtual scalpel displacement x(t) at a constant velocity. As such, this model represents an elliptical excision process as a virtual hybrid system with K linear regimes, in which the blade velocity \(\dot{x}(t)={v}_{k}\) is feedback-regulated by means of switching between the discrete regimes v1,…,vK. In this work, we show that such formulation can bring a greater insight into the analysis of surgical skill when compared to the descriptive statistics approach more commonly applied in this area. In the next section, we focus on the problem of inferring the parameters of our model from force measurements.

Excision as a switching linear dynamical system

The switching linear dynamical system30,31,32,33,34,35 is an example of a broader class of hybrid system, in which globally nonlinear dynamics are approximated by a series of linear systems. In the generative model of a switching linear dynamical system, the switching between each of its K linear regimes is described by a discrete hidden state variable st ∈ {1,…,K}. The evolution of st is characterized by K×K transition matrix Q that captures the probabilities of state transitions, i.e. P(st∣st−1). The continuous hidden state vector \({{{{{{{{\bf{z}}}}}}}}}_{t}\in {{\mathbb{R}}}^{D}\) evolves according to a D × D dynamics matrix A, and the observation vector \({{{{{{{{\bf{y}}}}}}}}}_{t}\in {{\mathbb{R}}}^{L}\) is generated according to an L × D observation matrix C, as follows:

$${{{{{{{{\bf{z}}}}}}}}}_{t}={{{{{{{{\bf{A}}}}}}}}}^{(k)}{{{{{{{{\bf{z}}}}}}}}}_{t-1}+{{{{{{{{\bf{w}}}}}}}}}_{t}^{(k)},$$

(4)

$${{{{{{{{\bf{y}}}}}}}}}_{t}={{{{{{{{\bf{C}}}}}}}}}^{(k)}{{{{{{{{\bf{z}}}}}}}}}_{t}+{{{{{{{{\bf{v}}}}}}}}}_{t}^{(k)}.$$

(5)

where A(k) and C(k) are associated with a regime st = k, and \({{{{{{{{\bf{w}}}}}}}}}_{t}^{(k)}\) and \({{{{{{{{\bf{v}}}}}}}}}_{t}^{(k)}\) are the disturbance and observation noise, respectively.

In this work, we model the elliptical excision process with two discrete linear regimes, k ∈ {L,U}. Each regime corresponds to a constant virtual velocity of the blade, and satisfies vL < vU (we call L — a lower regime, and U — an upper regime). For each of these linear regimes, we model the uncertainty in the constant velocity as \({\tilde{v}}_{k} \sim {{{{{{{\mathcal{N}}}}}}}}\left({v}_{k},{\sigma }_{k}^{2}\right)\), where \({\sigma }_{k}^{2}\) is the variance of the velocity noise in the regime k. The continuous hidden state vector \({{{{{{{{\bf{z}}}}}}}}}_{t}=\left[\begin{array}{c}{g}_{t}\\ {x}_{t}\\ 1\end{array}\right]\), comprises gt and xt, the latent cutting force and virtual displacement of the blade at time step t, respectively. Since we only measure the cutting force, the observable yt is a scalar that represents the force measurement at time step t. The continuous dynamics in the linear regime k is \({{{{{{{{\bf{A}}}}}}}}}^{(k)}=\left[\begin{array}{ccc}\alpha &\beta &0\\ 0&0&{\tilde{v}}_{k}\\ 0&0&0\end{array}\right]\), where constants α and β define the displacement-to-force relationship of the Maxwell model, and are found by transforming the transfer function, equation (3), into the equivalent state space form. The observation matrix in the linear regime k is \({{{{{{{{\bf{C}}}}}}}}}^{(k)}=\left[\begin{array}{ccc}\gamma &\delta &0\end{array}\right]\), where γ and δ are the observation constants from the state space representation of the Maxwell model’s transfer function. In this work, we set the spring constant E = 1 N cm−1 and the damping coefficient η = 0.5 N s cm−1, which yield α = − 2, β = 1, γ = − 2 and δ = 1 constant values. The parameters were chosen such that estimated displacements approximately match the actual distance travelled by the scalpel. Finally, given the uncertainty captured in the velocity \({\tilde{v}}_{k}\), we can further assume the disturbance-free dynamics (\({{{{{{{{\bf{w}}}}}}}}}_{t}^{(k)}\) is zero vector) and noise-free observations (\({{{{{{{{\rm{v}}}}}}}}}_{t}^{(k)}=0\)).

A graphical representation of this generative model is shown in Fig. 7a. There are several ways to infer the parameters of this class of models from observations. For example, the variational approach to learning in switching linear dynamical systems34 approximates the posterior probabilities of the hidden states by optimizing evidence lower bound. In this study, we bypass the inference of discrete hidden states st by assuming that velocities \(\dot{x}(t)\) are fully observable under the assumption of the Maxwell model (Fig. 7b). This turns the switching linear dynamical system inference into a problem of learning an HMM36, fully characterized by transition probability matrix Q (Fig. 7c) and the emission probabilities defined by vk and \({\sigma }_{k}^{2}\), for each of the linear regimes k. Given the virtual velocity profiles \(\dot{x}(t)\), this model can be easily fit using the Expectation-Maximization algorithm37.

a A graphical model representation of the generative model, where s is a discrete state, \(\dot{x}\) is blade’s velocity, x is blade’s displacement, g is excision force, y is force measurement and t is a time step. Shaded nodes represent the observed variables. b Hidden Markov Model (HMM) with a hidden discrete state st (the cutting regime at time step t), and an observable virtual velocity \({\dot{x}}_{t}\). c Markov chain with two cutting regimes defined by the transition matrix Q. d Model fitting (1, 2 and 3) and data generation processes (4, 5 and 6). (1) Actual measurements of forces collected during the trials. (2) The virtual displacement derived from the force measurements using the Maxwell model. (3) The virtual velocity profiles (finite differences of the displacement profiles) are used to train the HMM. (4) The velocity sampled from the trained HMM (blue line). (5) and (6) The synthetic displacement and force (blue lines), generated by the model.

Figure 7d provides an overview of the model fitting process. First, the virtual displacement profiles are derived from the force measurements using the inverse of the transfer function, specified by Maxwell model parameters, equation (3). Then, the obtained displacement profiles are numerically differentiated for estimation of the virtual velocities \(\dot{x}(t)\). Finally, the obtained virtual velocity profiles are used to fit an HMM with the Expectation-Maximization algorithm. (Examples of incision forces generated by the model when fit to each of the medical students are shown in Supplementary Fig. 8).

Sensorized cutting instrument

We constructed a uniaxial force sensor based on Texas Instrument’s LDC1612 inductance-to-digital converter (LDC) and a 3D printed flexible element. The LDC provides reliable position measurements at submicron resolution38, which in combination with a flexible element with a known stress-strain characteristic, enables the construction of displacement-based force sensors. The LDC measures the distance between a conductive target and an inductive coil using the resonant sensing principle. The inductive coil in parallel with the capacitor forms a resonant circuit in which the alternating current flowing through the inductor generates an alternating magnetic field. As a result of Faraday’s law, the alternating magnetic field induces eddy currents on the surface of the conductive target as a function of the target displacement. As per Lenz’s law, these eddy currents create an opposing magnetic field that reduces the nominal inductance of the resonant circuit, and hence, increases the resonant frequency. The LDC measures this frequency shift and thus provides information about the target’s displacement with respect to the inductor. By fixing the target to the free end of the flexure with a known stress-strain characteristic, a displacement measurement can be transformed into a force measurement.

The designed cutting tool consists of two key components, 1) a printed circuit board with an inductive coil, and 2) a flexure with a conductive target. The schematic for the uniaxial force sensor is shown in Supplementary Fig. 9b. The inductor is implemented as a circular planar coil of 8 mm diameter as shown in Supplementary Fig. 9a. In the rest configuration of the flexure, the effective 8.6 μH inductor (in parallel with 330 pF capacitor) focuses the alternating magnetic of 2.985 MHz frequency into the conductive target located 3.4 mm below. In our design, we used 10 mm square aluminium film of 0.2 mm thickness. The displacement range of the target is restricted to 1.6 mm, with a minimum distance to the inductor of 1.8 mm. When the flexure is at its maximum displacement configuration, the resonant frequency shifts from 2.985 MHz to 3.025 MHz (40 kHz shift, 1.3% of the nominal resonance at zero displacement). According to ref. 39, the maximum effective resolution achievable with the given frequency variation is 14-15 bits. The dimensions of the printed circuit board are 100 mm x 13.5 mm. The 4-layer board incorporates differential sensor coils, the LDC1612 inductance-to-digital converter, an MSP430F5528 microcontroller, power supply circuitry and a USB connector. The microcontroller configures the LDC via the I2C interface, implements USB Communication Device Class, processes and streams sensor data to a host computer.

The displacement is established by a one-piece 3D printed flexure, in which the free end displaces the conductive target under the presence of external force. As with any displacement-based force sensor, one of the main challenges is to maximize the stiffness of the flexure, while achieving the desired sensitivity. 3D printing provides a relatively easy way of experimenting with various design parameters, such as stiffness, strength, and geometry, as well as printing process parameters, such as material, printing orientation, etc. In this study, we use a blade flexure with design parameters shown in Supplementary Fig. 9b. The flexure was 3D printed with an Ultimaker 3 Extended printer using PLA thermoplastic, 0.2 mm layer height, 20% infill (triangle pattern) and 0.4 mm nozzle diameter. The extruder temperature was set to 205 °C, the travel speed was set to 70 mm per second and the perimeter layers were set to 3. The printing was done at room temperature controlled in a range between 19 and 21 °C. With these settings, the printed element was approximately 50 microns wider in XY direction.

Supplementary Fig. 9c shows the results of the incremental load test. During the test, a fully assembled device was incrementally loaded by ten 100 g calibrated weights (i.e. from 0.98 N to 9.8 N). The load was applied at the midpoint of the blade interface. The hysteresis (defined as the maximum difference between loading and unloading samples relative to the full-scale output) is 3.9%. The dotted line on the graph represents the linear least squares fit to the loading curve. The maximum deviation from the linear fit (non-linearity) is 1.4% of the full-scale output and the sensitivity of the sensor is 3752 counts per newton. Finally, the measured accuracy (maximum standard deviation of sensor output at the maximum measured load and relative to the maximum measured load, i.e. to 9.8 N) is 0.58%.

Tissue phantom

Supplementary Fig. 9d illustrates the design and material composition of the multilayered phantom used in this study. The design consists of a gelatin base that simulates the recoil of subcutaneous tissues, and a stack of three silicone layers that mimic the mechanical properties of human skin. The outer silicon layer is reinforced by pre-tensioned power mesh fabric that increases the tear strength of a sample. The gelatin base and silicon layers are coupled through a thin layer of an ultrasound gel. The fully assembled phantom has dimensions of 160 mm x 160 mm x 30 mm.

The fabrication of each phantom comprised of the following procedure. 64 g of gelatin powder (240 Bloom) was spread across 640 ml of cold water and left unstirred for 20 min, then simmered and stirred until fully dissolved. The liquid was poured into a 3D printed mould (160 mm x 160 mm x 25 mm volume container) wrapped in cellophane film and was left to solidify overnight in a refrigerator.

Next, a square piece of power mesh fabric (180 mm x 180 mm) was secured to the working surface under a slight amount of tension. 20 ml of two-part silicone rubber (Smooth-On EcoflexTM 00-30, shore hardness 30) was thoroughly mixed in a 1:1 ratio for 2 min and poured onto the center of stretched fabric in the series of three pours. The silicone-saturated mesh was then left for 45 min to cure. When cured, the next layer of 20 ml silicone (Smooth-On EcoflexTM GEL with shore hardness 000-35) was mixed and poured over. Finally, the second batch of 25 ml Smooth-On EcoflexTM 00-30 was poured over the pre-cured silicone layers. The silicone sample was left to cure for 4 h.

The cured silicone sample was placed on the full set gelatin base with a power mesh-reinforced layer presenting the skin surface. The remaining edges of the power mesh are trimmed to match the surface area of the phantom. The fully assembled phantom is stored in a refrigerator prior to each experiment.

The design of the phantom was selected after extensive validation with a single experienced surgeon, and selected for its realistic viscoelastic properties. A total of seven phantom designs were evaluated according to the perceived realism of pressing, stretching, pinching and cutting the phantom surface. All evaluated designs consisted of a gelatin base with 100 g per litre concentration and varying combinations of silicone layers. We have chosen Smooth-On EcoflexTM Gel, Smooth-On EcoflexTM 00-30 and Smooth-On Dragon SkinTM (shore hardness 10A) silicone rubbers to represent very soft, soft and hard phantom layers, respectively. Supplementary Table 5 shows the phantom design ranking (from least to most realistic). A few summary points:

Softer silicone rubbers (shore hardness < 30) appear more realistic.

The combination of silicone layers with varying hardness increases realism. Single-layer designs were scored lowest, while three-layer designs were rated as most realistic.

The hardness gradient (with a harder outer layer) plays a role in the realism of shear loads (e.g. stretching the skin).

The hardness of the bottom layer plays role in pressing load and can mimic the age of the skin.

Statistics and reproducibility

The statistical analysis was performed using open-source Python libraries SciPy (https://scipy.org/) and Pingouin (https://pingouin-stats.org/build/html/index.html). The elliptical excision force model was trained using open-source Python package hhmlearn (https://hmmlearn.readthedocs.io/en/stable/). For reproducibility, all data processing, analysis, modeling and figure generation routines were written using Jupyter Notebook.

[ad_2]

**This article was originally published by a ****www.nature.com ****. Read the ** **Original article here.**