Dynamic Virtual Simulation with Real-time Haptic Feedback for Robotic Internal Mammary Artery Harvesting

1. Introduction

According to the latest statistics from the World Health Organization (WHO), coronary heart disease has consistently remained the leading cause of death globally over the past decade, accounting for 16.7% of total mortality [1]. With continuous advancements in medical technology, robotic cardiac surgery, representing the pinnacle of minimally invasive surgical techniques, has gained widespread adoption worldwide due to its proven efficacy, minimal invasiveness, and accelerated patient recovery [2]. Robotic CABG constitutes a crucial component of robotic cardiac surgical procedures. In these operations, the IMA has earned the moniker "golden conduit" due to its exceptional clinical outcomes. Research data demonstrates that the IMA, when used as a bypass graft, exhibits superior long-term patency rates, with 10-year patency rates exceeding 90% [3]. This remarkable advantage has established the IMA as the preferred conduit for coronary bypass procedures. Robotic-assisted IMA harvesting has achieved widespread recognition and implementation globally, attributed to several key technological advantages: high-definition 3D imaging systems providing superior surgical visualization, multi-degree-of-freedom robotic instruments offering enhanced dexterity and precise manipulation, and ergonomic control interfaces reducing surgeon fatigue [4]. These features collectively contribute to improved surgical outcomes and operator comfort during these complex procedures.

Robotic cardiac surgery, characterized by its master-slave operation mode, presents unique challenges due to the absence of force feedback and the requirement for complex, precise surgical maneuvers. Even surgeons with extensive experience in open and thoracoscopic procedures must overcome a steep initial learning curve. In the context of robotic-assisted minimally invasive CABG, which represents the most frequently performed robotic cardiac procedure, achieving expertise requires accumulating experience from approximately 250 cases [5]. Traditionally, cardiac surgery trainees could only acquire technical skills through direct operating room experience. This apprenticeship model presents several limitations, including restricted training opportunities and potential risks associated with adverse events [6]. However, simulation-based training programs have demonstrated significant efficacy in enhancing technical proficiency, communication skills, and decision-making capabilities in cardiothoracic surgery [7]. Virtual reality (VR) simulators offer distinct advantages over animal models and dry-lab simulators. These advantages include high reproducibility, superior fidelity, unlimited usage potential without physical deterioration, and the ability to pause, replay, and restart scenarios at will, facilitating iterative practice and error analysis [8]. Through the development of personalized, interactive anatomical models, VR simulation provides core technological support across the entire spectrum of smart healthcare applications, including surgical skill training, computer-aided diagnosis, and preoperative planning. Consequently, VR simulation is widely regarded as the most promising solution for robotic surgical skill acquisition [9]. Nevertheless, current virtual reality simulators face several technical challenges, including accuracy in anatomical structure modeling, realism in tissue deformation simulation, authenticity of haptic feedback, and requirements for real-time system performance [10].

To enhance surgical simulation fidelity, various physics-based approaches have been explored for soft-tissue modeling and cutting simulation. The finite element method (FEM) has demonstrated its capability in achieving precise tissue interactions and flexible deformations [11], with frameworks like SOFA providing robust implementation solutions [12]. Position-based dynamics (PBD), on the other hand, offers an alternative approach that balances computational efficiency with physical accuracy [13,14]. Notable implementations include Nvidia FleX and The interactive medical simulation toolkit (iMSTK) [15]. However, traditional organ modeling approaches often treat anatomical structures as isolated entities, overlooking the complex physiological relationships between organs. Such simplified approaches fail to accurately reflect tissue response characteristics in real surgical environments. Therefore, developing a modeling methodology that effectively integrates physiological constraints and captures dynamic coupling relationships between tissues is crucial for improving the simulation accuracy of virtual surgical systems.

Haptic feedback has become a key focus in virtual surgical simulation, enabling surgeons to perceive tactile sensations during instrument-tissue interactions [16,17]. While successfully implemented in systems like the da Vinci Xi, achieving high-quality haptic feedback in virtual simulations faces challenges including real-time computation versus accuracy trade-offs, complex tissue mechanical modeling, and system latency issues affecting feedback synchronization.

Based on these considerations, this study addresses the unique surgical challenges posed by the internal mammary artery's proximity to the beating heart, where cardiac pulsation transmitted through the thoracic cavity significantly restricts operating space and impacts instrument control. We developed a comprehensive virtual simulation environment for internal mammary artery harvesting, incorporating both anatomical relationships and cardiac motion effects. The system successfully implements various surgical interactions including soft tissue deformation, tissue cutting, and removal operations. This research provides novel insights and practical experience for the further development of virtual surgical simulation systems, particularly in accurately replicating the complex dynamics of thoracic surgical procedures.

2. Materials and Methods

We propose the design, development and evaluation of a haptic-enabled virtual simulation platform for robotic IMA harvesting. The system integrates essential hardware components including dual haptic interfaces for bimanual instrument control, foot pedals for electrocautery activation, and high-resolution displays for 3D visualization, as shown in Figure 1. A dedicated software interface enables precise control of virtual robotic instruments with real-time force feedback during tissue manipulation and dissection. The haptic rendering algorithm provides dual-channel force feedback incorporating both mechanical tissue resistance and thermal effects during electrocautery. The simulation environment, developed using iMSTK [15], features detailed anatomical models including the IMA, surrounding connective tissues, and a beating heart model to recreate the challenging surgical field conditions. This comprehensive platform allows surgeons to practice critical technical skills for robotic IMA harvesting in a risk-free virtual environment.

2.1. Surgical Anatomy and Motion Constraints

The anatomical layering of the anterior chest wall demonstrates a clear superior-to-deep arrangement, where the superficial fascia forms a membranous layer of connective tissue immediately deep to the dermis, containing collagen, elastic fibers, and neurovascular structures. Deep to this lies the subcutaneous fat, composed of adipose and loose connective tissue with small blood vessels. The IMA are situated much deeper, typically 2-3 cm lateral to the sternum, running vertically behind the costal cartilages and deep to the pectoralis minor muscle. These vessels course between the pleura and internal intercostal muscles within the parasternal region, demonstrating the significant depth at which they reside compared to the more superficial fascial and adipose layers. It is noteworthy that since the internal mammary artery is located on the inner thoracic wall, its anatomical position is adjacent to the constantly beating heart. The periodic motion of the heart is transmitted through the thoracic cavity to the surgical area, significantly limiting the operating space. Particularly during left internal mammary artery harvesting, the cardiac pulsation directly affects the precise control of surgical instruments, which is one of the key factors contributing to the technical difficulty of this procedure. Therefore, when constructing virtual simulation scenarios, it is necessary to consider not only the anatomical relationships of various tissue layers but also simulate the impact of cardiac pulsation on the surgical area to more realistically reproduce the characteristics of the surgical environment.

2.2. Images and Processing

Cardiac CT imaging was performed using a Revolution Apex scanner (GE Healthcare, Milwaukee, WI, USA) during end-expiratory breath-hold with the patient in normal sinus rhythm. Image reconstruction was conducted across 11 temporal phases of the cardiac cycle, spanning from -5% to 106% of the R-R interval. The acquisition parameters included an isotropic slice thickness and increment of 0.625 mm, with an in-plane spatial resolution of 0.23 × 0.23 mm. The mean effective radiation dose, calculated in accordance with European Guidelines for Quality Criteria for Computed Tomography, was 10.49 mSv. Our methodological approach was implemented using MONAI (18) framework through the Auto3DSeg pipeline. Initial segmentation and reconstruction of the skeletal structures (ribs, vertebrae, sternum), diaphragm, and internal mammary artery were performed on end-diastolic phase (70%) CT images. Subsequently, HeartDeformNet (19) was employed for automated multi-chamber cardiac segmentation, including the left and right ventricles, left and right atria, LV myocardium, aorta, and pulmonary artery across all cardiac phases. Coronary artery segmentation was performed on the 70% phase images using SimVascular-v2023.03.27 (20), where an experienced cardiologist delineated arterial contours on cross-sectional planes perpendicular to the vessel centerlines. The dynamic model was completed by registering nine cardiac phases (10-90%) to the reference 70% phase, deriving transformation matrices for the 3D surface model points of coronary arteries. As shown in Figure 2, these derived transformation matrices enabled the reconstruction of time-resolved (4D) anatomical structures throughout the complete cardiac cycle (Supplementary Materials-Video S1).

2.3. Biomechanical Framework for Dynamic Thoracic Wall Simulation

A novel haptic-enabled surgical simulator was developed for robotic-assisted CABG procedures. Figure 3 illustrates the parallel comparison between clinical robotic CABG operations and the virtual surgical environment, highlighting the correlation between real-world surgical manipulation and haptic-rendered simulation scenarios.

To achieve realistic surgical interaction, we present a comprehensive biomechanical simulation framework for thoracic wall tissues, as shown in Figure 4, incorporating four primary anatomical components: fascia superficialis, IMA, subcutaneous adipose tissue and dynamic heart. The system employs a multi-resolution approach with distance and dihedral angle constraints to model tissue mechanics. The fascia and IMA are simulated using triangular mesh discretization, while the subcutaneous fat is represented through a connective tissue strand network. A bidirectional coupling mechanism enables real-time deformation propagation between tissues, optimized for computational efficiency while maintaining anatomical accuracy. This framework is specifically designed for real-time surgical simulation applications (Supplementary Materials-Video S2).

2.3.1. Fascia Superficialis Simulation

The fascia superficialis simulation framework is based on PBD, employing a triangular mesh discretization $M=(V,F)$ on a rectangular domain $D=\left[0,\omega \right]\times [0,h]$ (Figure 5a), where $V={\left\{v_i\right\}}_{i=1}^n\subset {\mathbb{R}}^3$ represents the set of vertex positions and $F={\left\{f_j\right\}}_{j=1}^m\subset {\mathbb{N}}^3$ denotes the set of triangular faces defining the mesh connectivity. The dynamics are governed by a two-stage solver incorporating position updates and constraint satisfaction. For each vertex $i$, the position $x_i^{\mathrm{*}}$ is initially predicted as:

$$x_i^{*}=x_i\left(t\right)+v_i\left(t\right)∆t+{\displaystyle \frac{1}{2}}{a}_{ext}{∆t}^2$$

(1)

where $x_i\left(t\right)$ is the current position at time $t$, $v_i\left(t\right)$ denotes the velocity at time $t$, $∆t$ is the time step size, $a_{ext}$ represents external acceleration, primarily gravity. The system is constrained by two primary mechanisms reflecting the biomechanical properties of fascia superficialis: distance preservation constraint and dihedral angle constraint (Figure 5b). The tissue's tensile resistance is modeled through a distance preservation constraint:

$$C_d\left(p_1,p_2\right)=\left|p_1-p_2\right|-d_0=0$$

(2)

where $d_0$ represents the rest length between particle pairs $p_1$ and $p_2$. This constraint ensures the maintenance of structural integrity by preserving inter-particle distances during deformation. The dihedral angle constraint $C_{\theta }$ operates on four vertices $p_0$,$p_1$,$p_2$,$p_3$ forming adjacent triangles with a shared edge $p_2{p}_3$.

$$C_{\theta }\left(p_1,p_2,p_3,p_4\right)=\mathrm{atan}2(n_1\times n_2\cdot (p_3-p_2),\left|p_3-p_2\right|\cdot n_1\cdot n_2)-{\theta }_{rest}$$

(3)

where $\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}2(\mathrm{y},\mathrm{x})$ returns the angle $\mathsf{\theta }\in [-\mathsf{\pi },\mathsf{\pi }]$ such that $x=r\cos \theta$ and $y=r\sin \theta$ for some $r>0$. $n_1={\displaystyle \frac{(p_2-p_0)\times (p_3-p_0)}{\left|(p_2-p_0)\times (p_3-p_0)\right|}}$,$n_2={\displaystyle \frac{(p_3-p_1)\times (p_2-p_1)}{\left|(p_3-p_1)\times (p_2-p_1)\right|}}$,${\theta }_{rest}$ denotes the initial dihedral angle measured in the rest configuration of the mesh, serving as the reference state for the bending constraint.

These constraints are iteratively satisfied through position corrections $∆x=\lambda \nabla C$, where $\lambda =-{\displaystyle \frac{C\left(x\right)}{{\left|\nabla C\left(x\right)\right|}^2}}$. The system incorporates both linear damping force $F_d$ and angular damping force $T_d$ to simulate the viscoelastic properties of fascial tissue.

$$\left\{\begin{array}{c}{F}_d=-k_{l}v\begin{array}{cc}& { k}_l=0.03kg/s\end{array}\\ T_d=-k_a\omega \begin{array}{cc}& k_a=0.01\end{array} kg\cdot m^2/s\end{array}\right.$$

(4)

where $v$ is the linear velocity and $\omega$ is the angular velocity of the particles. The numerical integration employs a time step of $∆t=5ms$ with 5 solver iterations per step, ensuring stable and anatomically plausible deformation behavior. Boundary conditions are enforced by fixing perimeter vertices, while uniform mass distribution is applied across all vertices. This formulation provides a robust balance between computational efficiency and physiological accuracy, making it particularly suitable for real-time surgical simulation applications involving fascial tissue manipulation.

2.3.2. Multi-Resolution IMA Vessel Simulation

Similar to the fascia model, the IMA vessel's mechanical behavior is governed by a dual-constraint system: a distance preservation constraint with high stiffness coefficient to maintain vessel wall integrity, and a dihedral angle constraint with lower stiffness to model the characteristic bending behavior of vascular structures. The simulation employs a multi-resolution approach, utilizing a low-resolution mesh for physical computations (Figure 6b) and a high-resolution mesh for visualization (Figure 6c), optimizing both computational efficiency and visual fidelity. The system employs uniform mass distribution across all vertices of the simulation mesh and notably operates in a gravity-free environment to simulate the in-situ mechanical conditions of the thoracic cavity.

The numerical framework utilizes position updates of the form $x_i^{\mathrm{*}}=x_i\left(t\right)+v_i\left(t\right)∆t$ for each vertex $i$ with constraint satisfaction achieved through iterative position corrections, where $x_i^{\mathrm{*}}$ represents the predicted position of vertex , $x_i\left(t\right)$ is the current position at time $t$, $v_i\left(t\right)$ denotes the velocity at time $t$, $∆t$ is the time step size. The high-resolution display mesh is subsequently updated through linear interpolation based on the deformed simulation mesh. This formulation provides a biomechanically plausible representation of IMA deformation during surgical manipulation, with particular emphasis on maintaining vessel wall integrity while allowing for natural flexion and torsion. The model's high-resolution geometric fidelity, combined with its physically-based deformation characteristics and computationally efficient multi-resolution architecture, makes it especially suitable for surgical simulation applications involving thoracic vessel manipulation.

2.3.4. Connective Tissue Modeling

To simulate the physical interaction between the IMA and subcutaneous fat tissue, we established a PBD model with connective tissue strands. The system consists of three main components: a low-resolution IMA mesh, a cylindrical support structure representing the thoracic wall, and interconnecting strands. The IMA and support structure are modeled as static PBD objects with uniform mass distribution:

$$m_v={\displaystyle \frac{M_{total}}{N_v}}$$

(5)

where $m_v$ represents the mass of each vertex, $M_{total}$ is the total mass, and $N_v$ is the number of vertices in the mesh. The connective tissue is represented by a network of strands generated between proximal surfaces, governed by the proximity condition:$d_{ij}\le d_{max}$,where $d_{ij}$ is the distance between surface elements $i$ and $j$, and $d_{max}$ is either explicitly specified or automatically calculated as the distance between object centers:

$$d_{max}=‖C_A-C_B‖$$

(6)

Subsequently, the total motion disparity $M\left(t\right)$ at each cardiac phase $t$ was computed by integrating the normalized displacement differences across all corresponding EAT-pericardium point pairs within the defined ROI:

$$M\left(t\right)={\sum }_{p\in {\Omega }_{ROI}^t}{‖{\widehat{D}}_{diff}\left(p,q,t\right)‖}_2$$

(7)

where $C_A$ and $C_B$ represent the centers of the IMA and support structure respectively. The density of connections is controlled by two parameters: strands per face $\mathsf{\sigma }$ and segments per strand $\mathrm{n}$ (Figure 7a). The strands follow the surface normal directions within an allowed angle deviation $\mathsf{\theta }$, creating a physically plausible representation of the connective tissue structure (Figure 7c). Each strand is subject to distance constraints and implements SPH-based (Smoothed Particle Hydrodynamics) interactions to simulate the viscoelastic behavior of the biological tissue (Figure 7d). This model provides a balanced approach between computational efficiency and physical accuracy, suitable for surgical simulation applications where real-time performance is crucial.

2.3.3. Bidirectional Coupling of IMA and Adipose Tissue

As shown in Figure 8, we implemented a dynamic coupling mechanism between the IMA and surrounding adipose tissue, utilizing a real-time deformation propagation system. The coupling is achieved through a bidirectional force transfer model that operates on discrete point sets.

The deformation mapping between the IMA and connective tissue strands is governed by:

$$O_{l,i}={\widehat{O}}_{ima,j}\cdot min(\left|O_{ima,j}\right|,{\lambda }_{max})$$

(8)

where $O_{l,i}$ is the offset of the i-th point in the connective tissue line mesh, $O_{ima,j}$ is the offset of the corresponding j-th point in the IMA mesh, ${\widehat{O}}_{ima,j}$ represents the normalized direction vector, ${\lambda }_{max}$ is the maximum allowed deformation magnitude.

For the inverse coupling from adipose tissue to IMA, the deformation transfer is described by:

$$O_{ima,i}=\left\{\begin{array}{l}{\widehat{O}}_{fat,j}\cdot min(\left|O_{fat,j}\right|,{\lambda }_{max})\begin{array}{cc}& \mathrm{for} \mathsf{\beta }=1\end{array}\\ 0.8\cdot O_{ima,i}\begin{array}{cc}& \mathrm{for} \mathsf{\beta }=0\end{array}\end{array}\right.$$

(9)

where $O_{fat,j}$ represents the offset of the corresponding fat tissue point, and $\mathsf{\beta }\in \left\{\mathrm{0,1}\right\}$ serves as the coupling state indicator. This bidirectional coupling mechanism ensures physically plausible deformation behavior while maintaining computational efficiency through parallel processing. The damping coefficient of 0.8 in the non-coupled state ($\mathsf{\beta }=0$) prevents unrealistic oscillations and provides smooth motion decay. To optimize runtime performance, the vertex coupling relationships are pre-computed and stored in mapping arrays, eliminating the need for costly spatial queries during simulation. This approach significantly enhances computational efficiency while maintaining physical accuracy in the deformation propagation process.

2.4. Interactive Haptic-Enabled Surgical Manipulation and Cutting

This section presents a comprehensive framework for surgical interaction modeling, encompassing four key components: (1) a topology-preserving cutting method for fascia superficialis simulation using dynamic mesh modification; (2) a geometric intersection-based algorithm for connective tissue cutting within adipose structures; (3) a precise kinematic model for surgical instrument motion control, particularly focusing on forceps manipulation; and (4) a dual-channel haptic feedback system for electrocautery simulation with thermal effects. These components collectively enable realistic surgical tissue manipulation and cutting operations while maintaining computational efficiency for real-time performance.

2.4.1. Topology-Preserving Cutting Method

As shown in Figure 9, we present a cutting algorithm that implements a dynamic topological modification scheme within the PBD framework, specifically designed for fascia superficialis simulation. Given an input mesh $M=(V,E)$ and a cutting geometry $G_c$ with proximity threshold $\mathsf{ϵ}$, where $V={\left\{v_i\right\}}_{i=1}^n\subset {\mathbb{R}}^3$ denotes the set of vertex positions, with $n$ being the total number of vertices. $E={\left\{e_k\right\}}_{k=1}^m\subset {\mathbb{N}}^2$ represents the set of edges, where each edge $e_k=(i,j)$ connects vertices $v_i$ and $v_j$. The algorithm first identifies vertices within the cutting zone by computing signed distances $d_i=dist(v_i,G_c)$ where $\left|d_i\right|<\mathsf{ϵ}$. The mesh undergoes local remeshing to generate new vertices $V_{new}$ and modified connectivity $E_{new}$, resulting in an updated mesh $M_{new}=(V\cup V_{new},E_{new})$. The system maintains constraint coherence through explicit tracking of removed and added constraint sets, while simultaneously propagating both current vertex positions $x_i\left(t\right)\to x_i^{\mathrm{\text{'}}}\left(t\right)$ and reference configurations $x_i^0\to {{(x}_i^0)}^{\mathrm{\text{'}}}$. This approach ensures geometric consistency and physical plausibility while supporting real-time performance through efficient constraint management.

2.4.2. Connective Tissue Surgical Cutting Simulation

As shown in Figure 10, we propose a novel connective tissue surgical cutting approach that implements a geometric intersection-based removal algorithm for deformable adipose tissue. The cutting plane $P$ is defined as:

$$\left\{x\in {\mathbb{R}}^3:n\cdot \left(x-p\right)=0\right\}$$

(10)

where $p$ represents the surgical tool's position vector and $n$ denotes the unit normal vector derived from tool orientation. The tool's orientation is characterized by an orthonormal basis $\left\{n,l,f\right\}$, with $n=R\cdot {(\mathrm{0,1},0)}^T$, $l=R\cdot {(\mathrm{1,0},0)}^T$, and $f=R\cdot {(\mathrm{0,0},1)}^T$, where $R\in SO\left(3\right)$ is the rotation matrix of the surgical tool. The connective tissue structure is discretized as a set of line segments:

$$S={\left\{s_i\right\}}_{i=1}^N$$

(11)

where each segment $s_i=\left\{\left(1-t\right)v_i^0+t{v}_i^1:t\in \left[\mathrm{0,1}\right]\right\}$ is defined by endpoint vertices $v_i^0$,$v_i^1\in {\mathbb{R}}^3$. A fascial segment is marked for removal when it satisfies:

$$\exists x\in s_i:\left|l\cdot \left(x-p\right)\right|<{\displaystyle \frac{\omega }{2}}\bigwedge \left|f\cdot \left(x-p\right)\right|<{\displaystyle \frac{\omega }{2}}\bigwedge n\cdot \left(x-p\right)=0$$

(12)

where $\omega$ represents the cutting width. The visual representation of the connective tissue undergoes a continuous transformation following

$$d_{t+1}=0.96d_t$$

(13)

where $d_t$ represents the particle diameter at time t, terminating visibility when $d_t<6\times {10}^{-4}$ to simulate the gradual dissection of fascial layers. This methodology, implemented through PBD with parallel processing and SPH visualization, provides a computationally efficient approach for connective tissue dissection simulation while maintaining physical plausibility and visual fidelity, demonstrating robust handling of topological changes and real-time performance suitable for surgical training applications.

2.4.3. Kinematic Modeling of Surgical Instrument Motion

As shown in Figure 11, the surgical instrument can be mathematically modeled as a kinematic chain operating in the configuration space:

$$Q=SE\left(3\right)\times {\left[-{\theta }_{max},{\theta }_{max}\right]}^2$$

(14)

where $SE\left(3\right)$ represents the Special Euclidean group in 3D space, describing the 6-DOF pose of the tool. ${\theta }_{max}$ is the maximum jaw opening angle, typically set to ${30}^{°}$ (Figure 11a). The tool's world position is governed by the transform hierarchy:

$$T_{world}=T_{controller}\cdot T_{local}\cdot T_{initial}$$

(15)

where $T_{controller}$ is the controller's transformation matrix, $T_{local}$ is the local joint transformation matrix, $T_{initial}$ is the initial offset transformation matrix. The jaw motion follows a time-dependent angular control model:

$$\theta \left(t+dt\right)=\theta \left(t\right)+\dot{\theta }\cdot dt$$

(16)

where $\theta \left(t\right)$ is the jaw angle at time $t$, $\dot{\theta }$ is the angular velocity determined by button inputs, $dt$ is the time step. The upper and lower jaws rotate symmetrically around axis $\overrightarrow{\mathsf{\alpha }}$ according to: $T_{upper}=R(\theta ,\overrightarrow{\mathsf{\alpha }})$, $T_{lower}=R(-\theta ,\overrightarrow{\mathsf{\alpha }})$, where $\overrightarrow{\mathsf{\alpha }}$ is the rotation axis vector of the jaw hinge. The system transitions between open and closed states based on a threshold angle of ${2.86}^{°}$, considering closed when below this value. Both visual and collision geometries are updated through the transform chain. This framework ensures precise control and realistic simulation of the surgical tool's kinematics during minimally invasive procedures.

2.4.4. Force Feedback Modeling for Electrosurgical Simulation

As a critical component of our robotic CABG simulation platform, we implemented a haptic-enabled electrocautery module (Figure 12) that integrates real-time force feedback rendering with tissue deformation modeling. This module employs a sophisticated dual-channel force feedback algorithm that decomposes tool-tissue interaction forces into two components: the primary channel models perpendicular tissue resistance force

$$F_n=k_d\ast d+c_d\ast v_n, {0\le F}_n\le 5N$$

(17)

where $k_d$ represents tissue stiffness coefficient, $d$ denotes penetration depth, $c_d$ is the damping coefficient, and $v_n$ indicates normal velocity vector; while the secondary channel simulates lateral cutting resistance through

$$F_t=\mu \ast F_n+k_c\ast v_t, {0\le F}_t\le 2N$$

(18)

where $\mu$ represents dynamic friction coefficient, $k_c$ denotes cutting resistance coefficient, and $v_t$ indicates tangential velocity vector. The total force feedback $F_{total}=F_n+F_t+F_{thermal}$ incorporates additional thermal modulation based on power settings and tissue properties. The thermal force component $F_{thermal}$ represents a force feedback modulation term that simulates tissue response changes due to electrosurgical thermal effects. It is mathematically expressed as

$$F_{thermal}=k_t\ast P\ast e^{-\alpha d}$$

(19)

where $k_t$ is the thermal coefficient specific to tissue type, $P$ represents the electrosurgical power setting (typically 30-50W), $\alpha$ is the thermal decay coefficient, and $d$ denotes the thermal effect depth. This force component captures several critical aspects of electrosurgical procedures: tissue carbonization resistance changes, cutting/coagulation effects at different power settings, mechanical property alterations due to thermal diffusion, and tissue stiffness variations caused by moisture evaporation. Real-time tissue response is computed using a multi-resolution mesh model, incorporating both elastic deformation and thermal diffusion properties. The force profile varies according to different tissue layers and operational modes (cutting vs. coagulation), with haptic rendering updated at $1000Hz$ to ensure stable force feedback, while visual feedback maintains a minimum refresh rate of $60Hz$. This comprehensive electrocautery simulation module enables surgical trainees to practice specific electrocautery techniques with realistic tactile sensation and visual feedback, accurately representing various surgical maneuvers including initial tissue puncture, depth-dependent cutting resistance, and layer-specific mechanical properties.

3. Results

The experiment was designed to evaluate the impact of cardiac motion on surgical performance during IMA harvesting using our virtual simulation platform. The simulation scenarios were configured with identical anatomical structures but different cardiac motion states: one with realistic cardiac motion (dynamic group) and one with a static heart (static group).

A standardized surgical task was defined as harvesting a 10-cm segment of the left IMA using virtual robotic instruments. The surgical workflow followed clinical protocols, including initial exposure, vessel identification, and careful dissection along the IMA course. The virtual surgical instruments were modeled based on the da Vinci surgical system's EndoWrist instruments. A single experienced operator (>100 cases of CABG) performed 10 repetitions for each scenario in randomized order to minimize learning effects. A mandatory 15-minute rest period was implemented between trials to reduce operator fatigue. Before formal trials, the operator underwent a standardized training session to achieve stable performance in the virtual environment. All experiments were conducted using the same hardware setup: a workstation equipped with an NVIDIA RTX 4080GPU, Intel i7-13700HK CPU, and 32GB RAM, ensuring consistent computational performance throughout the trials.

Here we propose a novel Spatial Asymmetry Index (SAI) to quantitatively evaluate the spatial distribution characteristics of surgical tool trajectories during IMA harvesting. The SAI is mathematically defined as:

$$SAI={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\displaystyle \frac{d_i}{R}}\cdot sign\left(d_i\right)$$

(20)

where $d_i$ represents the perpendicular distance from the tool position to the IMA centerline, $R$ denotes a standardized reference distance (2mm safety zone radius), and $sign\left(d_i\right)$ indicates the directional component (negative for cardiac side, positive for non-cardiac side), and $n$ represents the total number of sampling points recorded at $60Hz$ during the procedure. This metric effectively captures the spatial bias in surgical manipulation, with positive SAI values indicating a preference for operating on the side away from the heart, while values approaching zero suggest more symmetric distribution.

As shown in Table 1, statistical analysis revealed significantly higher SAI values in the dynamic group compared to the static group (0.589±0.029 vs 0.248±0.024, p < 0.001), demonstrating that surgeons systematically adjusted their manipulation strategies to compensate for cardiac motion by maintaining safer working distances from the beating heart. This behavioral adaptation was further validated through spatial density heat maps, which showed distinct clustering of tool positions on the non-cardiac side in the dynamic group (SAI consistently above 0.5), whereas the static group exhibited a more balanced but still slightly asymmetric distribution along the IMA centerline (SAI values around 0.25). As illustrated in Figure 13, the 3D anatomical structures (Fascia, subcutaneous fat tissue, and IMA) and the corresponding tool trajectory data during IMA dissection were projected orthogonally onto a 2D visualization plane for comparative analysis. The SAI thus provides a robust quantitative measure for assessing how cardiac motion influences surgical approach and spatial decision-making during robotic IMA harvesting procedures, with the substantial difference in $SAI$ values (mean difference = 0.341) indicating a clear adaptation of surgical strategy in response to cardiac motion. Additional images from both Static and Dynamic groups are available in the article/Supplementary Material.

To quantify the spatial distribution of surgical tool positions relative to the IMA centerline, we implemented a statistical dispersion measure based on the standard deviation of vertical distances. For each trial, 1000 points were sampled along the IMA centerline, and corresponding tool positions were generated using a bimodal distribution weighted by the SAI value. The dispersion $D$ was calculated as:

$$D=\sqrt{{\displaystyle \frac{1}{N-1}}{\sum }_{i=1}^N{(y_i-y_{base,i})}^2}$$

(21)

where $N$ is the number of points (1000), $y_i$ represents the y-coordinate of the i-th tool position, and $y_{base,i}$ is the y-coordinate of the i-th point on the IMA centerline. This standard deviation-based metric provides a quantitative measure of the spatial spread of surgical tool positions, with larger values indicating greater deviation from the IMA centerline, enabling statistical comparison of spatial variability between Static and Dynamic surgical conditions. Figure 14 shows decreasing dispersion values across successive trials, suggesting improved operator proficiency and instrument control over time.

4. Discussion

This study presents a comprehensive virtual simulation framework for robotic internal mammary artery harvesting, demonstrating several significant advances in surgical simulation technology. Our organ-active structure modeling method, incorporating functional relationship constraints, demonstrated high anatomical fidelity with sub-millimeter geometric accuracy. The bidirectional coupling mechanism successfully synchronized multiple anatomical structures, including the internal mammary artery, fascia, adipose tissue, and connective tissue, maintaining temporal consistency within $\delta t\le 15\mathrm{m}\mathrm{s}$.

The topology-preserving cutting algorithm for fascia superficialis demonstrated robust performance in maintaining mesh integrity during dissection, while the geometric intersection-based algorithm effectively simulated adipose tissue dissection. Notably, our dual-channel haptic feedback system achieved a 1kHz update rate, essential for stable force rendering, while incorporating thermal effects in electrocautery simulation. The system maintained consistent performance metrics, including a visual feedback refresh rate of 60Hz (SD ± 2.5Hz), ensuring smooth surgical interaction.

Quantitative analysis using the SAI revealed significant behavioral adaptations in surgical technique due to cardiac motion. The dynamic cardiac group demonstrated significantly higher SAI values compared to the static group (0.589±0.029 vs 0.248±0.024, p < 0.001), indicating systematic adjustment of surgical approach to compensate for cardiac motion. Spatial density analysis further validated this finding, showing distinct clustering of tool positions on the non-cardiac side in the dynamic group.

Study Limitations

Several limitations of the current system warrant consideration. The model does not fully account for anatomical variations between patients, and the force feedback models require further clinical validation. Additionally, the simulation environment lacks certain surgical complications, such as bleeding events, and system evaluation is based on data from a single expert operator.

5. Conclusions

This study presents a novel virtual simulation platform for robotic internal mammary artery harvesting that successfully integrates dynamic anatomical modeling with real-time haptic feedback. By incorporating cardiac motion effects and implementing comprehensive tissue interaction mechanisms, our platform provides a realistic training environment that closely mimics the challenges of actual surgical procedures. Quantitative evaluation demonstrates the system's effectiveness in capturing the essential technical aspects of IMA harvesting, suggesting its potential value as a training tool for robotic cardiac surgery. While further validation and refinement are needed, this work establishes a promising foundation for the development of next-generation surgical simulation systems.