Design and Modelling of a Continuum Robot for Distal Lung Sampling in Mechanically Ventilated Patients in Critical Care

1. Introduction

Critically ill patients who develop respiratory failure and require mechanical ventilation (MV) suffer a high morbidity and mortality. Indeed, Covid-19 patients who require MV, have a mortality approaching 40% in some case series. Once MV, patients are at high risk of developing secondary infections and other secondary complications. Rapid and accurate sampling of the distal lung is an important diagnostic procedure to guide therapeutic interventions. However, despite the high cost and attributable morbidity and mortality, diagnosis of diseases in the distal lung in MV patients is not standardised, lacks reproducibility and requires expert operators. Often, this leads to empirical treatments such as broad spectrum antibiotics which are then very difficult to deescalate, thus compounding the exposure of patients to non-targeted antimicrobials and promoting antimicrobial resistance. Pulmonary infiltrates in MV critically ill patients are a common occurrence and a major diagnostic challenge. Endobronchial secretions such as mucus and often hinder manually steered bronchoscopes, leading to poor sampling results. Hence, the aim of this paper is to develop a robotic bronchoscope that democratises sampling of the lung in MV ICU patients and enables non-skilled operators to safely sample disparate regions of the human lung to improve diagnostic accuracy and therapeutic interventions.

Bronchoscopy is a common diagnostic modality for the early detection of lung diseases (see Figure 1). During bronchoscopy, a thin tube (bronchoscope) is passed through the vocal cords into the airways to reach potential regions of the lung for directed sampling. Due to relatively large dimensions of the bronchoscope used for sampling (> 5 mm), bronchoscopy of MV patients is challenging. Another major drawback of the current technology is reliance on manual insertion, which is difficult due to the limited Degrees of Freedom (DoFs) of the bronchoscope, i.e., rotation and insertion.

FIGURE 1

Figure 1. A schematic of lung bronchoscopy in ICU, showcasing the insertion of the robotic bronchoscope through the mechanical ventilator and inside the lung.

To address the aforementioned challenges, we have developed a miniaturised continuum robot for lung bronchoscopy. A continuum robot has a continuously elastic structure and can traverse tightly curved 3D paths in confined spaces and reach desired positions deep inside human cavities. Continuum robots retain force transmission capability and offer great dexterity, thus, enabling optimal therapies when deeply seated pathologies are targeted (Burgner-Kahrs et al., 2015). Continuum robots have been explored for various interventions including laparoscopy (Wu et al., 2019), cardiac surgery (Fagogenis et al., 2019), neuro-surgery (Mattei et al., 2014), and eye surgery (Mitros et al., 2020).

The proposed bronchoscope is a continuum robot comprised of several parallel rods that can be bent via pushing/pulling of the rods. A continuum robot composed of several constrained push/pull rods is commonly known as a multi-backbone robot, first introduced in Gravagne and Walker (2000). Simaan et al. introduced the first surgical multi-backbone robot for dexterous tool manipulation in robotics surgery (Simaan et al., 2004; Ding et al., 2013). Xu et al. (2015) improved this design using a “dual continuum” actuation mechanism that increases modularity. Several researchers have explored the possibility of using a parallel multi-backbone approach without constraints, allowing more dexterous robots with increased DIFs per section (Bryson and Rucker, 2014; Wang et al., 2019). Multi-backbone robots have been commonly proposed for abdominal surgeries (Garbin et al., 2019; Riojas et al., 2019; Wu et al., 2019).

A major challenge in deployment of miniaturised continuum robots is real-time and precise modelling. There are several different kinematic and dynamic models presented in the literature (see Webster and Jones, 2010; Burgner-Kahrs et al., 2015 for a detailed review). The most common model for multi-backbone robots is a geometric model proposed in Simaan et al. (2004). The model has been used to control the motion of the robot as well as contact forces at the robot's tip (Bajo and Simaan, 2016). The geometric model assumes the robot curvature is constant and provides an accurate description of the robot's differential kinematics for large scale movements. However, due to the effects of unknown boundary conditions and the constant curvature assumption, the model's prediction of the robot shape and micro-scale movements are not accurate. To overcome this challenge, Del Giudice et al. (2017) proposed a method to improve micro-scale motion of a multi-backbone robot using modulation of the flexural rigidity of the rods. Another commonly method for modelling of multi-backbone robots is Cosserat rod theory (Bryson and Rucker, 2014; Wang et al., 2019). However, the Cosserat rod theory results in a relatively large boundary value problem (BVP) that should be solved for every rod in the robot and are computationally expensive. As a result, less accurate modelling methods are still attractive due to their low computational cost (Kaouk et al., 2014; Bajo and Simaan, 2016).

In this paper, we develop a bronchoscope using a miniaturised multi-backbone robot. The bronchoscope is mounted on a linear stage that can be used to automatically insert/retract the bronchoscope to reach targeted positions in the distal lung. Next, we develop a geometrically exact model of the robot that considers both the geometry of robot and mechanical properties of the backbones. The model results in a reduced order BVP and can be used to predict the shape of the bronchoscope without the constant curvature assumption. Furthermore, we develop a novel nonlinear observer that significantly improves the computational efficiency of the model to estimate the solution of the proposed model in real-time. Finally, simulations and experiments are performed to validate the design and the modelling Theory.

In the next section, section 2.1, the robot architecture and bronchoscope design is presented. Section 2.2 details the model of the bronchoscope. Section 2.4 outlines the detail of the observer design. In section 3, simulations and experimental results are performed to evaluate the design and quantify the accuracy and computational efficiency of the model. Concluding remarks appear in section 4.

2. Materials and Methods

2.1. System Design and Prototyping

This section describes the design and engineering of the robotic bronchoscope. The mechanical system design begins with DoFs discussion. To improve the dexterity of the bronchoscope, we propose a novel design that allows the robotic bronchoscope to bend in 3D at two points. The tip of the bronchoscope is composed of two segments shown in Figure 2. Each segment is actuated by 3 nitinol (NiTi) rods with an outer diameter of 0.475 mm which are passed through holes located on fixtures surrounding the bronchoscope. The circular fixtures are employed to avoid buckling of the rods. An additional silicone rod shown in blue in Figure 2 is acting as the main backbone. It has an outer and inner diameters of 2.3 mm and 2 mm, respectively and is rigidly connected to the fixtures to ensure they cannot move relative to each other. The fixtures' outer and inner diameters are 4.5 and 2.4 mm, respectively. Length of proximal segment at the tip of the bronchoscope is 40 mm, length of the distal segment is 500 mm, and the overall length of the bronchoscope is 540 mm. The end-effector is actuated via the push-pull of the 6 rods. In contrast to the cable driven bronchoscopes, the proposed design employs in-compressible Nitinol rods to offer more bending curvature via pushing of the rods. Furthermore, a 7th DoFs is employed for the insertion and retraction of the end-effector into the airways.

FIGURE 2

Figure 2. The robotic bronchoscope. (A) The inlet shows the tip of the bronchoscope which is composed of two segments that can be independently bent. By pulling/pushing the wires at each segment the bronchoscope can bend in 3D space. (B) The bronchoscope prototype placed inside a 3D printed lung model. An electromagnetic tracker (Aurora electromagnetic tracking system, NDI, Canada) is placed at the tip of the bronchoscope to measure its tip position in real-time. (C) Camera view from the endoscopic camera placed inside the working channel of the bronchoscope.

All DoFs are actuated by brushless DC motors (Maxon Motors) with a gearhead with a 150:1 reduction and a quadratic encoder. Each motor is controlled via a position controller module with a built in PID controller (EPOS4 Compact 50/5 CAN). The controllers employ the encoders feedback to accurately control the position of the motor shaft. The position controllers communicate with a PC via the CAN protocol. A CAN-to-USB interface (Kvaser Inc., CA, USA) is used to connect the position controllers to the PC.

The motors are connected to lead screws that convert the power generated by the motor into feed velocity for pulling/pushing the rods. The lead screws are carrying a v-shaped 3D printed part that is connected to the rods (shown in Figure 2) and travels along the lead screw to pull/push the rods. Additionally, 6 linear potentiometers are used to accurately measure the displacement of the rods.

Figure 2 shows the developed robot and an inlet showcases the different segments that the manipulator comprises.

2.2. Geometrically Exact Model of the Robot

We use the Cosserat-rod theory (Antman, 2005; Rucker et al., 2010) to model the robot. First, we present the model for a robot with one bendable segment. Next, we generalise the model to a robot with more segments. The following notation is used throughout the paper: x, x, and X denote a scalar, a vector, and a matrix, respectively. A complete summary of variables and operators is given in the Appendix. The symbols used are summarised in a nomenclature section.

A schematic of the robot is shown in Figure 3. The robot comprises a main backbone (shown in blue) rigidly connected to the fixtures and three NiTi rods (shown in red) fixed at the end fixture. The three rods can pass through the rest of the fixtures and have enough clearance to not create forces and moments but rather follow the curvature of the main backbone. The relative position of each rod with respect to the main backbone (d_i, i = 1, 2, 3 in Figure 3) is given by

$\begin{array}{ll}{d}_i={[\delta \text{cos}({\beta }_i),\delta \text{sin}({\beta }_i),0]}^T,& (1)\end{array}$

where δ is the rods' distance from the robots centroid (see Figure 3) and β_i is the relative angular position of each rod with respect to the main backbone

$\begin{array}{ll}{\beta }_i=\alpha +(i-1)\frac{2\pi }{3},i=1,2,3,& (2)\end{array}$

with α shown in Figure 3.

FIGURE 3

Figure 3. A schematic of the continuum robot with one bent segment. The main backbone is modelled as a Cosserat rod under external point force (F) and distributed load (f). The cross section view shows the position of the rods with respect to the main backbone.

The robot main backbone is modelled as a long, slender, one-dimensional Cosserat rod endowed with a Darboux frame attached to every point on its arc with the z axis of the frame tangent to the curve. The rod is under an external point force [F(t)] and distributed constant load (f) simulating the weight of the fixtures. The configuration of the rod can be defined using a unique set of 3D centroids, r(s, t) :[0, ℓ] × [0, ∞] → ℝ³ × [0, ∞], and a family of orthogonal transformations, R(s, t) :[0, ℓ] × [0, ∞] → so(3) × [0, ∞]. The position of the main backbone is defined by

$\begin{array}{ll}r\prime (s,t)=\text{R}(s,t){\mathit{e}}_3,{\text{R}}^{\prime }(s,t)=\text{R}(s,t){[\text{u}(s,t)]}_{\times },& (3)\end{array}$

where $u(s,t)={\left[u_x(s,t),u_y(s,t),u_z(s,t)\right]}^T$ is the curvature vector of the deformed backbone, [.]_× operator is the isomorphism between a vector in ℝ³ and its skew-symmetric cross product matrix, and $e_3={\left[0,0,1\right]}^T$ is the unit vector aligned with the z-axis of the global coordinate frame. Assuming the rods are made of linear elastic isotropic materials, we can derive the constitutive equations for calculating the instantaneous curvature of the rod (Rucker et al., 2010)

$\begin{array}{l}{\mathit{u}}^{\prime }(s,t)=-K^{-1}[{[\text{u}(s,t)]}_{\times }K\mathit{u}(s,t)+\\ {[{\mathit{e}}_3]}_{\times }{\text{R}}^T(s,t)(\mathit{F}(t)+(l-s)f)],& (4)\end{array}$

where l is the length of the main backbone and K = diag(EI, EI, GJ) is the stiffness matrix for the whole robot; E is the robot's Young's modulus; I is the second moment of inertia; G is the shear modulus; J is the polar moment of inertia. It is assumed that the cross section of the robot is symmetric and the products of inertia are negligible (i.e., I_xy = I_xz = I_yz ≃ 0)

In practice, the robot curvature u(s, t) and position r(s, t) are unknown and should be estimated as the function of the length of the three rods (ℓ_i, i = 1, 2, 3). We can estimate each rod's total arc length as

$\begin{array}{ll}{\ell }_i(t)={\displaystyle {\int }_0^l\Vert }{\mathit{r}}^{\prime }i(s,t)\Vert \text{d}s,& (5)\end{array}$

where ‖⋅‖ denotes the ℓ₂-norm and r_i(s, t) is the position of ith rod given by

$\begin{array}{ll}{r}_i(s,t)=r(s,t)+\text{R}(s,t)d_i.& (6)\end{array}$

Substituting (6) in (5) and simplifying the equations using (3) yields

$\begin{array}{ll}{\ell }_i(t)={\displaystyle {\int }_0^l\Vert {\mathit{e}}_3+[\text{u}(}s,t){]}_{\times }{\mathit{d}}_i\Vert \text{d}s,& (7)\end{array}$

Now, we can write the system of differential equations governing the motion of the robot using (3), (4), and (7)

$\begin{array}{ll}{\mathit{r}}^{\prime }(s,t)=\text{R}(s,t){\mathit{e}}_3,& (\mathrm{8a})\end{array}$

$\begin{array}{ll}{\text{R}}^{\prime }(s,t)=\text{R}(s,t){[\text{u}(s,t)]}_{\times },& (\mathrm{8b})\end{array}$

$\begin{array}{l}{\mathit{u}}^{\prime }(s,t)=-K^{-1}[{[\text{u}(s,t)]}_{\times }K\mathit{u}(s,t)+\hfill \\ {[{\mathit{e}}_3]}_{\times }{\text{R}}^T(s,t)(\mathit{F}(t)+(l-s)f)],\hfill \end{array}(\mathrm{8c})$

$\begin{array}{ll}{\ell }^{\prime }i(s,t)=\Vert {\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i\Vert ,i=1,2,3,& (\mathrm{8d})\end{array}$

with the following boundary conditions

$\begin{array}{ll}\mathit{r}(0,t)={\left[000\right]}^T,& (\mathrm{9a})\end{array}$

$\begin{array}{ll}\text{R}(0,t)=\text{I},& (\mathrm{9b})\end{array}$

$\begin{array}{ll}{u}_z(0,t)=0,& (\mathrm{9c})\end{array}$

$\begin{array}{ll}{\ell }_i(0,t)=0,& (\mathrm{9d})\end{array}$

$\begin{array}{ll}{\ell }_i(l,t)=L_i(t),i=1,2.& (9e)\end{array}$

The model defined by (8) and (9) accepts the overall length of the first two rods L_i, i = 1, 2 as inputs and predicts the robot curvature u(s, t), position r(s, t), and length of the third rod ℓ₃(t). We note that the length of the third rod is always defined by the length of the first and second rod.

Additionally, (8) and (9) form a boundary value problem. In the absence of external torques, the initial curvature of the robot along z direction is zero (9c). However, the initial curvatures along x and y directions [i.e., u_x(0, t) and u_y(0, t)] are unknown. In addition, the first and second rods' arc length ℓ_i(s, t), i = 1, 2 are defined both at the base (s = 0) and the tip of the robot (s = l) by (9d, 9e).

Moreover, the model given in (8) is quasi-static and solved using the separation of variables. To solve the equations, it is assumed that at a given time, time-dependent variables are constant and the equations are solved in spatial domain (with respect to s) using standard methods such as the Runge-Kutta or Adams-Bashforth families of algorithms. Shooting methods can be used to solve the boundary value problem. A shooting method consists of using a nonlinear root-finding algorithm to iteratively converge on values for u_x(0, t) and u_y(0, t), in order to satisfy (9e). Next, the time-dependent variables are updated [i.e., L_i(t)], and the equations are solved again in the spatial domain.

2.3. Multi-Segment Robot

Here, we generalise the model given in (8) for a multi-backbone robot with multiple bending segments shown in Figure 4. It is assumed that the robot is composed of n segments with lengths of l_j, j = 1, .., n. Each segment is actuated via 3 parallel rods fixed at the end the segment. Thus, there are n rods and the jth segment contains 3 × (n + 1 − j) rods.

FIGURE 4

Figure 4. A schematic of multi-backbone robot with multiple bending segments, dashed lines denote break points.

To model the robot, we start from the 1st segment containing n × 3 rods and use (8) to estimate the curvature, position of the main backbone, and the lengths of the cables up to the next segment. Next, at the junction where the segment ends (shown with dashed lines in Figure 4) we enforce the appropriate boundary conditions. The boundary conditions to be enforced across each transition point between sections are as follows: (1) The position and orientation of each tube must be continuous across the boundary, i.e.,

$\begin{array}{ll}\mathit{r}(s^{-},t)=\mathit{r}(s^{+},t),\text{ R}(s^{-},t)=\text{R}(s^{+},t),& (10)\end{array}$

(2) considering the static equilibrium and the fact that the rods apply no torque around z direction:

$\begin{array}{ll}{u}_z(s^{-},t)=u_z(s^{+},t),& (11)\end{array}$

(3) at the distal end of each segment, we have a boundary condition for the length of the rods that end:

$\begin{array}{ll}{\ell }_j(s,t)=L_j(t).& (12)\end{array}$

We repeat this process for the rest of the segments. We note that the curvatures along x and y at each break point are unknown. A shooting method must be used to iteratively converge on values for {u_x(0, t), u_y(0, t), u_x(l₁, t), u_y(l₁, t),…, $u_x({\sum }_{j=1}^n{l}_j,t),u_y({\sum }_{j=1}^n{l}_j,t)\}$, in order to acquire the desired length for the rods. Solving the BVP numerically is computationally intensive. The computational cost of the model is a significant obstacle in deployment of such designs and more efficient numerical methods are needed. To this end, we study the design of a novel observer that can rapidly estimate the solution of robot's model without the need to solve the BVP.

2.4. Rapid Solution of the Robot's Model

Our main goal in this section is to design an observer that employs measurements of ℓ_i(l, t) through time to estimate correct value of u_y(0, t) and u_x(0, t) and ensures the boundary conditions in (9) are satisfied without the need to solve the BVP iteratively. First, we transform the model in (8) into an observable form that simplifies the design of the observer. Next, we design an observer rule that guarantees exponential convergence of the solution of the observable model to the solution of the robot model given in by (8) and (9).

We define the vector of missing initial values as

$\begin{array}{ll}\check{u}(0,t)={[u_y(0,t),u_x(0,t)]}^T.& (13)\end{array}$

To realise the effect of the missing initial value [i.e., $\check{u}(0,t)$] on the solution of the equations, we define four auxiliary variables, namely,

$\begin{array}{ll}{{A}}_i(s,t):=\frac{\partial {\ell }_i(s,t)}{\partial \check{u}(0,t)},i=1,2,3& (\mathrm{14a})\end{array}$

$\begin{array}{ll}{B}(s,t):=\frac{\partial u(s,t)}{\partial \check{u}(0,t)},& (\mathrm{14b})\end{array}$

$\begin{array}{ll}{C}(s,t):=\frac{\partial \left[{\text{R}}^T(s,t)(\mathit{F}(t)+(l-s)f)\right]}{\partial \check{\mathit{u}}(0,t)},& (\mathrm{14c})\end{array}$

$\begin{array}{ll}{D}(s,t):=\frac{\partial \left({\text{R}}^T(s,t)\mathit{f}\right)}{\partial \check{\mathit{u}}(0,t)}.& (\mathrm{14d})\end{array}$

Using the new variables defined by (14), one can derive the following generalised model of the multi-backbone robot (see the Appendix for derivation)

$\begin{array}{ll}{\mathit{r}}^{\prime }(s,t)=\text{R}(s,t){\mathit{e}}_3,& (\mathrm{15a})\end{array}$

$\begin{array}{ll}{\text{R}}^{\prime }(s,t)=\text{R}(s,t){[\text{u}(s,t)]}_{\times },& (\mathrm{15b})\end{array}$

$\begin{array}{l}{\mathit{u}}^{\prime }(s,t)=-K^{-1}[{[\text{u}(s,t)]}_{\times }K\mathit{u}(s,t)+\\ {[{\mathit{e}}_3]}_{\times }{\text{R}}^T(s,t)(\mathit{F}(t)+(l-s)f)],& (\mathrm{15c})\end{array}$

$\begin{array}{ll}{{\ell }^{\prime }}_i(s,t)=\Vert {\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i\Vert ,i=1,2,3,& (\mathrm{15d})\end{array}$

$\begin{array}{ll}{{{A}}^{\prime }}_i(s,t)=\frac{-{({\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i)}^T{[{\mathit{d}}_i]}_{\times }}{\Vert {\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i\Vert }ℬ(s,t),i=1,2,3,& (15e)\end{array}$

$\begin{array}{l}{{B}}^{\prime }(s,t)=K^{-1}[{[K\mathit{u}(s,t)]}_{\times }{B}(s,t)-\\ {[\mathit{u}(s,t)]}_{\times }K{B}(s,t)-{[{\mathit{e}}_3]}_{\times }{C}(s,t)],& (\mathrm{15f})\end{array}$

$\begin{array}{l}{{C}}^{\prime }(s,t)=[{\text{R}}^T(s,t)(\mathit{F}(t)+(l-s)f){]}_{\times }{B}(s,t)-\\ [\mathit{u}(s,t){]}_{\times }{C}(s,t)-{D}(s,t),& (\mathrm{15g})\end{array}$

$\begin{array}{ll}{{D}}^{\prime }(s,t)={[{\text{R}}^T(s,t)\mathit{f}]}_{\times }{B}(s,t)-{[\mathit{u}(s,t)]}_{\times }{D}(s,t).& (\mathrm{15h})\end{array}$

Now, we provide a set of initial conditions for (15) that ensures the solution of the observer model in (15) exponentially converges to the solution of the boundary value problem defined by (8) and (9).

$\begin{array}{ll}\mathit{r}(0,t)={\left[000\right]}^T,& (\mathrm{16a})\end{array}$

$\begin{array}{ll}\text{R}(0,t)=\text{I},& (\mathrm{16b})\end{array}$

$\begin{array}{ll}{u}_z(0,t)=0,& (\mathrm{16c})\end{array}$

$\begin{array}{ll}\left[\begin{array}{c}{u}_x(0,t)\\ u_y(0,t)\end{array}\right]=-{\displaystyle {\int }_0^t}{\left[\begin{array}{c}{{A}}_1^T(l,t)\\ {{A}}_2^T(l,t)\end{array}\right]}^{†}\text{P}\left[\begin{array}{c}{\ell }_1(l,t)-L_1(t)\\ {\ell }_2(l,t)-L_2(t)\end{array}\right]\text{ d}t,& (\mathrm{16d})\end{array}$

$\begin{array}{ll}{\ell }_i(0,t)=0,& (16e)\end{array}$

$\begin{array}{ll}{{A}}_i(0,t)=0,i=1,2,& (\mathrm{16f})\end{array}$

$\begin{array}{ll}{B}(0,t)=\left[100;010\right],& (\mathrm{16g})\end{array}$

$\begin{array}{ll}{C}(0,t)=0,& (\mathrm{16h})\end{array}$

$\begin{array}{ll}{D}(0,t)=0,& (\mathrm{16i})\end{array}$

where P is a symmetric positive definite matrix and ^† denotes the pseudo-inverse operator. We note that (16f-16i) are calculated based on the definition of the auxiliary variables in (14). (16d) is the main observer rule that guarantees the convergence of the observer (see the Appendix).

The observer given in (15) is quasi-static, similar to the robot's model in (8). However, instead of using an iterative BVP solver, it can be solved as an initial value problem using the initial values given in (16). At a given time t, time-dependent variables are assumed constant and the equations are solved in spatial domain. Next, the time-dependent variables are updated [i.e., L_i(t), u_x(0, t), u_y(0, t)]. The updated time-dependant variables are used to solve the equations in the spatial domain again. The observer can be generalised to a multi-segment robot following the approach discussed in section 2.3.

In the next section, series of simulation and experiments are performed to evaluate the model's accuracy and demonstrate the computational efficiency of the observer in comparison with common BVP numerical solvers.

3. Results

Simulations and experiments are performed to evaluate the proposed design and modelling theory. The bronchoscopic robot used in the simulations and experiments consists of two bendable segments, shown in Figure 2. The length of the first segment is 500 mm, and the length of the second segment (at the tip) is 40 mm. The outer diameter of the robot is 4.5 mm and the inner diameter of the robot is 2 mm. Twenty-seven circular fixtures each weighting 5 g were equally spaced along the length of the bronchoscope and were rigidly fixed to the main backbone shown in blue in Figure 2.

We performed experiments to identify the developed model parameters and validate the model. First we performed experiments to identify the model parameters. For the identification phase, each rod was commanded to either push or pull the end disks by 5 mm, making the robot to randomly bend to 12 different positions. We estimated the 3D shape of the robot using calibrated stereo rig comprising two Logitech HD Pro C922 webcams. The cameras were running at 1, 080p resolution. As identified through calibration using on average 30 views of a checkerboard, a single pixel corresponded to 0.25 × 0.25 mm on the image plane. Following calibration, the entry point of the robot, i.e., $s=0$ was estimated in 3D space via triangulation. The robot coordinate frame was aligned to a planar calibration target always visible by the cameras during the experiments.

Furthermore, manual backbone segmentation established the base and shape of the bronchoscope relative to the aligned calibration grid. Matching backbone points were selected in both images, and then triangulated to provide the 3D point cloud. This process is shown in Figure 5. The mean error of the 3D triangulation algorithm was equal to 6 pixels corresponding to 1.5 mm. The extracted 3D backbones were used to identify for the robot model parameters, namely, Young's modulus, E and shear moduli G of the robot and initial displacement of the rods, δℓ_i, i = 1, …, 6. The parameters were identified by fitting the kinematic model given in (8) to the shape of the robot estimated via the cameras at 12 different configurations. The identified parameters of the model and the known parameters of the model are given in Table 1.

FIGURE 5

Figure 5. Estimating robot's backbone shape using two calibrated cameras.

TABLE 1

Table 1. Physical parameters of the robot.

In the next step, to validate the model accuracy we commanded the robot to move to 20 different positions. The shape of the robot was estimated using the calibrated cameras and was compared to the shape of the robot predicted by the identified model. Figure 6 shows representative results. Results of the measurements including maximum, mean, and standard deviation of error of the model in predicting the robot tip position and the root-mean-squared error of the model in predicting robot shape are listed in Table 2. The root-mean-squared error is calculated as

$\begin{array}{ll}\text{RMSE}=\sqrt{\frac{{\displaystyle {\sum }_{j=1}^m}{\left(\Vert \widehat{r}-r{\Vert }_j\right)}^2}{m}},& (17)\end{array}$

and is used as a measure of the differences between the actual shape of the robot, $\widehat{r}$, and the model predicted shape, i.e., r, for m = 30 data points along the robot backbone.

FIGURE 6

Figure 6. A comparison of experimental bronchoscope's shape with model prediction at four different configurations.

TABLE 2

Table 2. Experimental results.

In the experiments, the robot tip was capable of bending up to 100° with respect to its backbone (see Figure 6). The maximum error of the model in estimating the position of the robot tip is 26.2 mm, corresponding to 1.9% of the robot's length.

Finally, we performed simulations to compare the computational efficiency of the observer with various shooting methods used to solve BVPs. Shooting methods consists of using a nonlinear root-finding algorithm to iteratively converge on values for u_x(0, t) and u_y(0, t) for each segment, in order to satisfy the boundary conditions (9), i.e.,

$\begin{array}{ll}\begin{array}{c}\text{Minimize}: Error:=\Vert \left[\begin{array}{c}{\ell }_1(l_1,t)-L_1(t)\\ {\ell }_2(l_1,t)-L_2(t)\\ {\ell }_4(l_1+l_2,t)-L_4(t)\\ {\ell }_5(l_1+l_2,t)-L_5(t)\end{array}\right]\Vert ,\\ \text{w}.\text{r}.\text{t}.: u_x(0,t),u_y(0,t),u_x(l_1,t),u_y(l_1,t).\end{array}\hfill & (18)\hfill \end{array}$

We compared the observer predictions with shooting method algorithms that employ three different root-finding algorithms, which to the best of authors knowledge, are the most commonly used BVP solvers. These solvers are:

1. Interior-point method (Byrd et al., 2000),

2. Quasi-Newton method with BFGS Hessian estimation (Curtis and Que, 2015),

3. Nelder-Mead method (Powell, 1973).

In the simulations, we pulled and pushed the cables from −5 to 5 mm at a frequency of 2π/10 Hz. The simulation runs for 10 seconds at sampling frequency of 50 Hz. The observer gain P used in the simulations was set to 70 × I, as this value was found to achieve the minimum prediction error. The optimally tolerance for all the algorithms were set to 10⁻⁶. The simulations are performed in Matlab on an Intel Core i7 (2.93 GHz) machine with 16 GB memory.

Figure 7 shows the robot's trajectory estimated via different methods. The observer and the Nelder-Mead method gave the best accuracy. The other two methods, namely, interior-point and quasi-Newton, gave substantial error at two points across the robot trajectory. Also, it can be seen that the observer has an error at the first sampling time but rapidly converges to the correct solution.

FIGURE 7

Figure 7. A comparison of bronchoscope's tip trajectory calculated by solving the robot's model using four different methods. The bronchoscope's backbone is shown at several configurations along the trajectory.

Figure 8 shows the error of the solvers and the observer in satisfying the boundary conditions given in (18). The observer error is the same order as the BVP solvers. The BVP solvers occasionally fail in minimizing the error, while the observer consistently maintains an error below 10⁻⁴ mm. Figure 8 compares the computational efficiency of the BVP solvers and the observer in terms of the time that each method takes to compute the solution of the model at each sampling time. Evidently, the observer is much faster than the BVP solvers and has lower standard deviation. The average time that the observer takes to estimate the model's solution is 0.0108 s, which is significantly faster than other solvers. We performed 10 more simulations, where, the robots rods are pulled/pushed at frequencies between π/5 Hz and π/50. The results are summarised in Table 3. The results demonstrate that the observer maintain similar error as the BVP solvers, while exhibiting superior computational efficiency. The average time that the observer takes to estimate the model's solution is 47 times faster than the fastest BVP solver, namely, the interior-point method.

FIGURE 8

Figure 8. A comparison of (A) accuracy and (B) computational efficiency of the observer with common BVP solvers. On each box in (B), the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points and the outliers are plotted individually using plus symbol.

TABLE 3

Table 3. Experimental results.

4. Discussion

In this paper, we presented the concept for and the design of a continuum robot for pulmonary endoscopy in MV patients. MV patients are at high risk of developing secondary infections and there is a need for a reliable and controlled sampling of the distal lung to guide therapeutic interventions. Current methods for diagnosis of diseases in the distal lung in MV patients are not standardised, lack reproducibility and require expert operators. Here, we proposed a novel robotic bronchoscope that can be used to democratise lung sampling and improve the accuracy and reliability of distal lung sampling in MV patients. The proposed design of the system considers the limitations and constraints of current bronchoscopy, i.e., limited dexterity, low repeatability, and relatively large size of the bronchoscope.

One of the main challenges in current bronchoscopy is navigating the tightly curved architecture of bronchial tree. Several studies (Coppola et al., 1998; Ulusoy et al., 2016) have reported that bifurcation angles of the bronchial tree including sub-carinal angles and inter-bronchial angles vary between 30 and 100°. The experimentally measured maximum bending angle of the proposed robotic bronchoscope is 100° with respect to the robot's main backbone, which enables the robot to traverse the tightly curved structure of airways. We note that the maximum bending angle and curvature of the robot is a function of the robot's interaction with the environment. In the future, we will study the bending capability of the robot in lung models to fully verify the effectiveness of the robot in navigating the bronchial tree.

The external diameter of traditional bronchoscopes is generally 5–6 mm with a working channel with inner diameter of 2 mm (Vachani and Sterman, 2008). The developed prototype is comparable with current technology and has an outer diameter of 4.5 mm with a working channel with inner diameter of 2 mm. Moreover, the bronchoscope is highly dexterous and possesses 7 DoFs in total. The continuum manipulator is able to bend in 3D at 2 different points along its backbone thanks to 6 push/ pull NiTi rods. The extra dexterity offered by the proposed design can potentially extend the reach of the clinical bronchoscopy.

One of the aims of this research is to democratise bronchoscopy in MV patients in the ICU via automating the procedure. To this end, we have proposed a new theoretical framework to model the robot that can be used in closed-loop control of the bronchoscope motion. Our novel mechanics-based model of the robotic bronchoscope can predict the shape of the robotic bronchoscope under external forces with an accuracy corresponding to 1.9% of its arc-length. We note that for long, slender continuum robots, tip error is highly dependent on the total arc length (Rucker et al., 2010) and robot's backbone's interaction with its surrounding environment. We note that this error can be further reduced via closed-loop control of the robot tip. A closed-loop controller can employ sensory feedback from the robot tip position to minimise the bronchoscope error in navigating the lung. In practice, electromagnetic trackers are placed at the tip of the bronchoscope to measure its tip position in real-time. The proposed design offers a 2 mm working channel that can be used to place such trackers, allowing real-time monitoring of robot position for closed-loop control.

Furthermore, we have demonstrated that our numerical framework can estimate the model's solution 47 times faster than the fastest existing solvers, enabling applications in real-time robotic control. Future work will focus on developing a closed-loop control strategy that uses the model and the feedback of the robot tip position acquired with electromagnetic trackers, to minimise the error of the robot tip in following a desired trajectory for sampling.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

MK conceived the study. ZM contributed to modelling and design of the robot. BT contributed to the experimental study. All authors contributed to manuscript writing, revision, and read and approved the submitted version.

Funding

This work was supported by the Medical Research Council [grants MR/T023252/1, and 8532390].

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Antman, S. (2005). Nonlinear Problems of Elasticity, 2nd Edn. Dordrecht: Springer.

Bajo, A., and Simaan, N. (2016). Hybrid motion/force control of multi-backbone continuum robots. Int. J. Robot. Res. 35, 422–434. doi: 10.1177/0278364915584806

CrossRef Full Text | Google Scholar

Bryson, C. E., and Rucker, D. C. (2014). “Toward parallel continuum manipulators,” in IEEE International Conference on Robotics and Automation (Hong Kong), 5264–5267. doi: 10.1109/ICRA.2014.6906943

CrossRef Full Text | Google Scholar

Burgner-Kahrs, J., Rucker, D. C., and Choset, H. (2015). Continuum robots for medical applications: a survey. IEEE Trans. Robot. 31, 1261–1280. doi: 10.1109/TRO.2015.2489500

CrossRef Full Text | Google Scholar

Byrd, R. H., Gilbert, J. C., and Nocedal, J. (2000). A trust region method based on interior point techniques for nonlinear programming. Math. Prog. 89, 149–185. doi: 10.1007/PL00011391

CrossRef Full Text | Google Scholar

Coppola, V., Vallone, G., Coscioni, E., Coppola, M., Maraziti, G., Alfinito, M., et al. (1998). Normal value of the tracheal bifurcation angle and correlation with left atrial volume. La Radiol. Med. 95, 461–465.

PubMed Abstract | Google Scholar

Curtis, F. E., and Que, X. (2015). A quasi-Newton algorithm for nonconvex, nonsmooth optimization with global convergence guarantees. Math. Prog. Comput. 7, 399–428. doi: 10.1007/s12532-015-0086-2

CrossRef Full Text | Google Scholar

Del Giudice, G., Wang, L., Shen, J., Joos, K., and Simaan, N. (2017). “Continuum robots for multi-scale motion: micro-scale motion through equilibrium modulation,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (Vancouver, BC), 2537–2542. doi: 10.1109/IROS.2017.8206074

CrossRef Full Text | Google Scholar

Ding, J., Goldman, R. E., Xu, K., Allen, P. K., Fowler, D. L., and Simaan, N. (2013). Design and coordination kinematics of an insertable robotic effectors platform for single-port access surgery. IEEE/ASME Trans. Mechatron. 18, 1612–1624. doi: 10.1109/TMECH.2012.2209671

PubMed Abstract | CrossRef Full Text | Google Scholar

Fagogenis, G, Mencattelli, M, Machaidze, Z, Rosa, B, Price, K, and Wu, F Autonomous robotic intracardiac catheter navigation using haptic vision. Sci. Robot. (2019) 4:eaaw1977. doi: 10.1126/scirobotics.aaw1977.

PubMed Abstract | CrossRef Full Text | Google Scholar

Garbin, N., Wang, L., Chandler, J. H., Obstein, K. L., Simaan, N., and Valdastri, P. (2019). Dual-continuum design approach for intuitive and low-cost upper gastrointestinal endoscopy. IEEE Trans. Biomed. Eng. 66, 1963–1974. doi: 10.1109/TBME.2018.2881717

PubMed Abstract | CrossRef Full Text | Google Scholar

Gravagne, I. A., and Walker, I. D. (2000). “Kinematic transformations for remotely-actuated planar continuum robots,” in Proceedings 2000 ICRA, Millennium Conference, IEEE International Conference on Robotics and Automation, Symposia Proceedings (Cat. No.00CH37065), Vol. 1 (San Francisco, CA), 19–26. doi: 10.1109/ROBOT.2000.844034

CrossRef Full Text | Google Scholar

Kaouk, J. H., Haber, G.-P., Autorino, R., Crouzet, S., Ouzzane, A., Flamand, V., et al. (2014). A novel robotic system for single-port urologic surgery: first clinical investigation. Eur. Urol. 66, 1033–1043. doi: 10.1016/j.eururo.2014.06.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattei, T. A., Rodriguez, A. H., Sambhara, D., and Mendel, E. (2014). Current state-of-the-art and future perspectives of robotic technology in neurosurgery. Neurosurg. Rev. 37, 357–366; discussion: 366. doi: 10.1007/s10143-014-0540-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitros, Z., Sadati, S., Seneci, C., Bloch, E., Leibrandt, K., Khadem, M., Cruz, L., et al. (2020). Optic nerve sheath fenestration with a multi-arm continuum robot. IEEE Robot. Autom. Lett. 5, 4874–4881. doi: 10.1109/LRA.2020.3005129

CrossRef Full Text | Google Scholar

Powell, M. J. D. (1973). On search directions for minimization algorithms. Math. Prog. 4, 193–201. doi: 10.1007/BF01584660

CrossRef Full Text | Google Scholar

Riojas, K. E., Anderson, P. L., Lathrop, R. A., Herrell, S. D., Rucker, D. C., and Webster, R. J. (2019). A hand-held non-robotic surgical tool with a wrist and an elbow. IEEE Trans. Biomed. Eng. 66, 3176–3184. doi: 10.1109/TBME.2019.2901751

PubMed Abstract | CrossRef Full Text | Google Scholar

Rucker, D. C., Jones, B. A., and Webster, R. J. III (2010). A geometrically exact model for externally loaded concentric-tube continuum robots. IEEE Trans. Robot. 26, 769–780. doi: 10.1109/TRO.2010.2062570

PubMed Abstract | CrossRef Full Text | Google Scholar

Simaan, N., Taylor, R., and Flint, P. (2004). “A dexterous system for laryngeal surgery,” in IEEE International Conference on Robotics and Automation, Vol. 1 (New Orleans, LA), 351–357. doi: 10.1109/ROBOT.2004.1307175

CrossRef Full Text | Google Scholar

Ulusoy, M., Uysal, I. I., Kıvrak, A. S., Ozbek, S., Karabulut, A. K., Paksoy, Y., et al. (2016). Age and gender related changes in bronchial tree: a morphometric study with multidedector CT. Eur. Rev. Med. Pharmacol. Sci. 20, 3351–3357.

PubMed Abstract | Google Scholar

Vachani, A., and Sterman, D. H. (2008). “Chapter 12 – bronchoscopy,” in Clinical Respiratory Medicine, 3rd Edn, eds R. K. Albert, S. G. Spiro, and J. R. Jett (Mosby, PA: Elsevier), 177–196. doi: 10.1016/B978-032304825-5.10012-1

CrossRef Full Text

Wang, J., Ha, J., and Dupont, P. E. (2019). “Steering a multi-armed robotic sheath using eccentric precurved tubes,” in 2019 International Conference on Robotics and Automation (ICRA) (Montreal, QC), 9834–9840. doi: 10.1109/ICRA.2019.8794245

PubMed Abstract | CrossRef Full Text | Google Scholar

Webster, R. J. III., and Jones, B. A. (2010). Design and kinematic modeling of constant curvature continuum robots: a review. Int. J. Robot. Res. 29, 1661–1683. doi: 10.1177/0278364910368147

CrossRef Full Text | Google Scholar

Wu, Z., Li, Q., Zhao, J., Gao, J., and Xu, K. (2019). Design of a modular continuum-articulated laparoscopic robotic tool with decoupled kinematics. IEEE Robot. Autom. Lett. 4, 3545–3552. doi: 10.1109/LRA.2019.2927929

CrossRef Full Text | Google Scholar

Xu, K., Zhao, J., and Fu, M. (2015). Development of the sjtu unfoldable robotic system (surs) for single port laparoscopy. IEEE/ASME Trans. Mechatron. 20, 2133–2145. doi: 10.1109/TMECH.2014.2364625

CrossRef Full Text | Google Scholar

A. Nomenclature

Appendix

Derivation of Auxiliary Variables

Here, we derive the differential equations governing the evolution of auxiliary variables given in (14), beginning with (14a).

Using (8d) and the chain rule of differentiation we have

$\begin{array}{l}{{A}}^{\prime }i(s,t)=\\ \frac{\partial }{\check{\mathit{u}}(0,t)}{\left[{({\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i)}^T({\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i)\right]}^{1/2}=\\ \frac{-{({\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i)}^T{[{\mathit{d}}_i]}_{\times }}{\Vert {\mathit{e}}_3+{[\text{u}(s,t)]}_{\times }{\mathit{d}}_i\Vert }{B}(s,t)\end{array}$

In deriving (15) we used the following identity

$\frac{\partial ({[\mathit{a}]}_{\times }\mathit{b})}{\partial \mathit{c}}=-{[\mathit{b}]}_{\times }\frac{\partial \mathit{a}}{\partial \mathit{c}}+{[\mathit{a}]}_{\times }\frac{\partial \mathit{b}}{\partial \mathit{c}}.$

We now use (8c) to derive the equations for calculating ${B}(s,t)$.

$\begin{array}{l}{{B}}^{\prime }(s,t)=\frac{\partial {\mathit{u}}^{\prime }(s,t)}{\check{\mathit{u}}(0,t)}={\text{K}}^{-1}[{[\text{K}\mathit{u}(s,t)]}_{\times }{B}(s,t)-\\ {[\mathit{u}(s,t)]}_{\times }\text{K}{B}(s,t)-{[{\mathit{e}}_3]}_{\times }{C}(s,t)]\end{array}$

Furthermore, ${C}(s,t)$ can be computed by taking the transpose of (8b), multiplying both sides by F(t) + (l − s)f, subtracting R^T(s, t)f, and finally taking partial derivative of both sides with respect to $\check{u}(0,t)$.

$\begin{array}{l}{{C}}^{\prime }(s,t)=[{\text{R}}^T(s,t)(\mathit{F}(t)+(l-s)f){]}_{\times }{B}(s,t)-\\ [\mathit{u}(s,t){]}_{\times }{C}(s,t)-{D}(s,t).\end{array}$

We can calculate ${D}(s,t)$ in a similar way. First, we take the transpose of (8b). Next, we multiply both sides by f. Finally, taking partial derivative of both sides with respect to $\check{u}(0,t)$ gives

$D^{\prime }(s,t)={[{\text{R}}^T(s,t)\mathit{f}]}_{\times }ℬ(s,t)-{[\mathit{u}(s,t)]}_{\times }D(s,t).$

Proof of Convergence and Stability

The error of the observer in satisfying the robot's model's boundary condition in (9) (i.e., robot's rods' lengths) is

$\begin{array}{l}ϵ(t)=\left[\begin{array}{c}{\ell }_1(l,t)\\ {\ell }_2(l,t)\end{array}\right]-\left[\begin{array}{c}{L}_1(t)\\ L_2(t).\end{array}\right]\end{array}$

To prove that this error exponentially converges to zero, we select the following Lyapunov candidate

$\begin{array}{ll}V=\frac{1}{2}{ϵ}^Tϵ.& (\mathrm{A1})\end{array}$

Taking the time derivative of V and replacing $\dot{ϵ}$ using (14a) we obtain

$\dot{V}={ϵ}^T\dot{ϵ}={ϵ}^T\left[\begin{array}{c}{A}_1^T(l,t)\\ A_2^T(l,t)\end{array}\right]\stackrel{.}{\stackrel{̌}{u}}(0,t).$

Substituting $\dot{\check{u}}(0,t)$ using (16d) gives

$\begin{array}{l}\dot{V}=-{\mathit{ϵ}}^T\left[\begin{array}{c}{{A}}_1^T(l,t)\\ {{A}}_2^T(l,t)\end{array}\right]{\left[\begin{array}{c}{{A}}_1^T(l,t)\\ {{A}}_2^T(l,t)\end{array}\right]}^{†}\text{P}\mathit{ϵ}\\ =-{\mathit{ϵ}}^T\text{P}\mathit{ϵ}\end{array}$

In deriving the above equation we used the following identity

$\begin{array}{l}\text{a}{\text{a}}^{†}=\text{a}{\text{a}}^T{(\text{a}{\text{a}}^T)}^{-1}=\text{I}.\end{array}$

P is symmetric positive definite. Thus, ∀t > 0, $\dot{V}$ is negative definite. Therefore, as t → ∞, $\dot{ϵ}(t)\to 0$. Consequently, the solution of the observer converges to the solution of the model given in (8).