The MUSHA hand is the main and innovative part of the teleoperation system, while the PSM of the da Vinci robot is already a known system. Then, in what follows we describe the multi-functionality of the MUSHA hand and its kinematic model which is crucial in the development of the entire control approach. According to Liu et al. (2019), and as anticipated in Section 1, the tool is a three-fingered underactuated miniature hand, each finger is composed of 12 segments and connected with a one DoF wrist, the kinematic model is composed of 37 joints connected by 36 links. An optoelectronic force sensor is designed, calibrated and presented in Liu et al. (2020), and located at the tip of each finger. In the range of the sensors, equal to 4 N, the maximal error on x-, y-, and z-direction is 0.28, 0.24, and 0.75 N, respectively. The MUSHA is actuated by the four motors of da Vinci PSM plus two additional motors. As a result, six joints are independent and directly connected to actuators by tendons, while the other ones are passive. For these reasons, it is possible to define a reduced kinematic model composed of actuated joints and equivalent links (Figure 2). The kinematic structure and the actuation system are designed to be multi-functional allowing the MUSHA hand performing three main methods of grasp, depicted in Figure 3: 1) power grasp (Po), 2) precision grasp (Pr) and 3) retractor (Rt). During the power grasp method, all fingertips and phalanges are involved to hold firmly tissue, in the precision grasp method only the fingertips are used to pinch tissue, and the retractor method allows the retraction of the tissues that come in contact with phalanges of index and middle fingertips. The relationship between joints and the tips position for a single robotic finger results as: x r , ρ = cos q 1 , ρ a 2 ⁡ cos q 2 , ρ + a 3 ⁡ cos q 2 , ρ + q 3 , ρ y r , ρ = sin q 1 , ρ a 2 ⁡ cos q 2 , ρ + a 3 ⁡ cos q 2 , ρ + q 3 , ρ z r , ρ = a 2 ⁡ sin q 2 , ρ + a 3 ⁡ sin q 2 , ρ + q 3 , ρ ∀ ρ ∈ T , I , M (1) where q _i,ρ is the value of the i-th joint of the robotic hand r ρ -th finger and a _i is the length of the i-th link. We refer with T, I, M the thumb, index, and middle finger, respectively. Thanks to the mechanical properties of the robotic tool tendon-driven transmission, q _3,ρ and q _2,ρ are decoupled and thus they can be independently controlled, it follows that the model can be rewritten as x r , ρ = cos q 1 , ρ a 2 ⁡ cos q 2 , ρ + a 3 ⁡ cos q 3 , ρ y r , ρ = sin q 1 , ρ a 2 ⁡ cos q 2 , ρ + a 3 ⁡ cos q 3 , ρ z r , ρ = a 2 ⁡ sin q 2 , ρ + a 3 ⁡ sin q 3 , ρ ∀ ρ ∈ T , I , M . (2)

The complete kinematic model of the MUSHA hand, p _r = K(⋅), is obtained exploiting (2) for each finger, thus p _r = K(q _1,T , q _2,T , q _3,T , q _1,I , q _2,I , q _3,I , q _1,M , q _2,M , q _3,M ). According to the sub-actuation constrains we define the wrist angle position, q _W , and the index and medium proximal one, qPIP,: q W = q 1 , I q W = q 1 , T + 2 π 3 q W = q 1 , M − 2 π 3  and  q P I P = q 2 , M q P I P = q 2 , I (3) we can simplify the kinematic relationship, obtaining: p r = K q W , q 2 , T , q 3 , T , q P I P , q 3 , I , q 3 , M (4) where p _r is a vector containing the fingertips position expressed in Cartesian coordinate. It is worth noticing that in the exploited model, the length of the links a ₂ and a ₃ depends on the joint values. Indeed, the coupling motion of the phalanx segments (generated by the tendons) causes links shrinkage. In accordance with the chord theorem (see Figure 4), the length of the links can be computed as a i , ρ = 2 r i , ρ ⁡ sin β i 2 . (5)

Let q i j , ρ be the joint j in the phalanx i of the finger ρ, ∀i ∈ {2, 3} and ∀ρ ∈ {T, I, M}. Consequently, it is possible to write the length of the chord in kinematic model and the length of the link in the actuated model as follows: a i , ρ = 2 L i , ρ ∑ j = 2 6 q i j , ρ + ∑ j = 1 5 q i j , ρ sin ∑ j = 2 6 q i j , ρ + ∑ j = 1 5 q i j , ρ 2 L i , ρ  if  q i , ρ = 0 (6) where L _i , ρ is the length of the arc of the i-th phalanx of the finger ρ. Assuming that the action of the tendons is equal on all the segments, the angular displacements are equal to a fraction of the equivalent joint in the reduced model: q i j , ρ = q i , ρ 6   ∀ j ∈ 1,6 ,  ∀ i ∈ 2,3  and  ∀ ρ ∈ T , I , M . (7)

As a consequence, ∀ρ ∈ {T, I, M} the single finger model becomes: x r , ρ = cos q 1 , ρ a 2 , ρ ⁡ cos q 2 , ρ + a 3 , ρ ⁡ cos q 3 , ρ (8) y r , ρ = sin q 1 , ρ a 2 , ρ ⁡ cos q 2 , ρ + a 3 , ρ ⁡ cos q 3 , ρ z r , ρ = a 2 , ρ ⁡ sin q 2 , ρ + a 3 , ρ ⁡ sin q 3 , ρ a 2 , ρ = 6 L 2 , ρ 5 q 2 , ρ sin 5 q 2 , ρ 6 L 2 , ρ  if  q 2 , ρ = 0 a 3 , ρ = 6 L 3 , ρ 5 q 3 , ρ sin 5 q 3 , ρ 6 L 3 , ρ  if  q 3 , ρ = 0

Therefore, the explicit full hand model (reported in the Supplementary Material) has 14 equations nine of which linearly independent, in six variables. Finally, we can express p ̇ r = J r w q ̇ r , where J r w is the Jacobian matrix of the robotic fingers of MUSHA expressed in reference frame of the wrist. The differential model is fully reported in the Supplementary (see Supplementary Section S5.2).

In this section, we report the simplified kinematic model of the hand, which will be instrumental in the design of the control technique. Even if a complete human hand model has about 30 DoF (Lee and Kunii, 1995), for the sake of simplicity and without loss of generality we use a simplified kinematic hand structure. We modeled each finger as a planar kinematic chain, with three 1-D hinges. As in Lisini et al. (2017), we assume that each finger has the metacarpal (MC) bone fixed with respect to the hand frame, and it is characterized by three DoFs. In a more complete model of the human hand, the thumb has at least five DoFs: two in the trapeziometacarpal (TM) joint, two in the metacarpophalangeal (MCP) joint, and one in the interphalangeal (IP) joint. However, the thumb abduction/adduction motion range of the MCP joint usually can be neglected so the thumb can be modeled with two DoFs. Index, middle, ring, and pinky fingers have one DoFs in the MCP revolute joint for flexion/extension, one in the proximal IP (PIP), and one in the distal IP (DIP). Figure 5 shows the model of the hand used in this work. Therefore, the exploited hand model has 20 DoFs. i.e., 14 DoFs describing the fingers and six DoFs (position and orientation) for the palm. Finally, we introduce the relationship between fingertips and joints velocities, which can be described by the differential model: p ̂ ̇ h = J h q h q ̇ h (9) where the matrix J h ∈ R 15 × 20 is the hand Jacobian, that maps the joint velocities ( q ̇ h ∈ R 20 ) to the Cartesian fingertip velocities ( p ̂ ̇ h ∈ R 15 ).

Starting from the definition of the human hand and the MUSHA hand models, the developed method to control the MUSHA hand and the da Vinci robot arm are explained. A standard point-to-point mapping technique is used to control the PSM, while an innovative approach based on the theory of mapping in object domain to teleoperate the MUSHA fingers is designed and implemented. We assume that the human and the robotic hands are constantly in touch with a similar virtual object, all possible contact points are called active points. The positions of MUSHA active points p _r are defined as the positions of fingertips, while the human active points ( p h ⊂ p ̂ h ) are given by the positions of palm (P), thumb, index, and middle fingertips: p h = p h , P , p h , T , p h , I , p h , M T . (10)

We indicate with p _h,ρ a vector of Cartesian positions of the ρ-th finger of the human hand h expressed in the world reference frame: p h , ρ = x h , ρ , y h , ρ , z h , ρ T , ∀ ρ ∈ P , T , I , M . (11)

Depending on the grasping method, we define the vector of interaction points, c _h , equal to a subset of active points of the human hand involved in the desired grasp, and evaluated as: c h = ∫ t c ̇ h τ d τ (12) being c ̇ h = I Ω p ̇ h  with  I ∈ R 12 × 12 (13) and Ω = 1 12 × 1 if method = Pr 0 1 × 3 , 1 1 × 9 T if method = P o 1 1 × 3 , 0 1 × 3 , 1 1 × 6 T if method = R t (14)

It is worth noticing that Ω is a selection vector that is evaluated on-line to allow switch control reference without discontinuity and according to the grasping method required. The palm-thumb distance is used to recognize the current grasping method.

The reference of PSM in cartesian space is calculated as integration over time τ of the velocity components of the user hand centroid c ̃ ̇ h until the time instant t. In this way, during a temporary disconnection between the patient side and control side, the relocation of the palm in more ergonomic poses or in a position far from the boundaries of the vision cone of the camera system is allowed. This results in: p P S M t = p P S M t − τ + ∫ t c ̃ ̇ h τ d τ (15) Where the value of centroid changes according to method of grasping, in fact: c ̃ h = ∑ ρ ∈ P , T , I , M c h , ρ N  with  N = 3  if method ≠ P o 4  if method = P o (16)

Similarly, the orientation of the human hand is exploited to orient the PSM. To this end Euler angles in roll (α), pitch (β), and yaw (γ) convention are computed as: α = 1 2 tan − 1 2 μ 0,1,1 μ 0,2,0 − μ 0,0,2 , (17) β = 1 2 tan − 1 2 μ 1,0,1 μ 0,0,2 − μ 2,0,0 , (18) γ = 1 2 tan − 1 2 μ 1,1,0 μ 2,0,0 − μ 0,2,0 , (19) being μ _i,j,k the second order central moments of cloud of contact points μ i , j , k = ∑ ρ ∈ P , I , T , M x h , ρ − x ̃ h j y h , ρ − y ̃ h i z h , ρ − z ̃ h k (21)

We exploited a sphere as virtual object to easily project the human hand motions into MUSHA fingers and wrist movements.

The roto-translations and deformations of a sphere can be completely defined by the following seven parameters: 1) o h ∈ R 3 Cartesian coordinate of the center of sphere; 2) θ h ∈ R 3 the orientation of the frame attached to the sphere with respect to the world frame, expressed as roll, pitch, and yaw angles; 3) r h ∈ R the sphere radius. Combining (13) and the aforementioned seven parameters, it is possible to obtain a linear relationship between human hand motions and sphere transformation c ̇ h , ρ = o ̇ h + θ ̇ h × c h , ρ − o h + r ̇ h c h , ρ − o h . (22)

Thus, there exists a Jacobian matrix J _hs such that its pseudo-inverse, J h s † , projects the human fingertip velocities into derivative parameters of the sphere: o h ̇ θ h ̇ r h ̇ = J h s † o h , c h I Ω p h ̇ (23)

The resulting Jacobian matrix depends on the current position of interaction points of human hand and on the sphere center that is a function of the active points of the human hand. Indeed, the sphere center is the circumcenter of a triangle with three of the current interaction points as vertices. o h = c h , i + | Δ h , k i | 2 Δ h , j i × Δ h , k i × Δ h , j i 2 | Δ h , j i × Δ h , k i | 2 + | Δ h , j i | 2 Δ h , k i × Δ h , j i × Δ h , k i 2 | Δ h , j i × Δ h , k i | 2 (24) where Δ _h,ki is the difference between the position vector of the k-th fingertip c _h,k in contact with the sphere and the i-th fingertip c _h,i , then, (23) can be rewritten as: o h ̇ θ h ̇ r h ̇ = J h s † p h I Ω p h ̇ (25)

As a further step towards the mapping, the sphere manipulated by the human is properly transformed into the one manipulated by the MUSHA through a scale factor K _s , that is equal to the ratio between the radius of human and MUSHA manipulated spheres. The algorithm computes the Jacobian of the sphere manipulated by the robot (i.e., J _rs ) and evaluates the velocities of the robotic fingertips, p ̇ r = J r s p r K s J h s † p h I Ω p h ̇ . (26)

Velocities are estimated by a closed loop kinematic inversion q ̇ r = J r w † q r J r s p r K s J h s † p h I Ω p h ̇ + K c p ̃ r − p r (27) where J r w , J _rs and J _hs are the Jacobian matrices of the robot, of the sphere manipulated by the robot, and of the sphere manipulated by the user, respectively. K c ( p ̃ r − p r ) is the corrective action on the error between the current position of MUSHA interaction points p _r and the desired ones p ̃ r , weighted by the gain K _c . It is worth pointing out that the pseudo-inversion of J _hs induces an indeterminacy in the estimated positions far from the current pose in case of multiple solutions, for this reason, the right position of the thumb in the retractor method may be not correctly computed. This situation is visually depicted in Figure 6. A projection into the null-space of mapping actions is used to overcome this limitation. To this end, (27) is extended as follows, q ̇ r = J r w † q r J r s p r K s J h s † p h I Ω p h ̇ + K c p ̃ r − p r + J r w I 6 × 6 − J r s K s J h s I Ω † J r s K s J h s I Ω u (28) where u = 0 6 × 1 if method ≠ R t K R t q R t − q r if method = R t (29) being K _Rt the gain that multiplies the error between current joint position q _r and q _Rt , the expected joint position in retractor pose. Regarding the initial condition, at the beginning of any task, the MUSHA lies in a rest position with all the fingers closed with joint configuration q r 0 = 0 6 × 1 . This particular configuration results in a kinematic singularity, in fact when q(t) = q _r0 the rank of Jacobian of robot decreases, R a n k ( J r w ) < 6 . For this reason, an initialization procedure is required. In accordance with q r 0 = arg max q det J r w q J r w T q q i , ρ ∈ q i , ρ min , q i , ρ max   ∀ i ,  ∀ ρ (30) we select q r 0 = [ 0 , π 9 , − π 9 , π 9 , − π 9 , − π 9 ] as initial position. It corresponds to the joint configuration with the maximum manipulability. In conclusion, combining (28–30), we obtain the entire human-MUSHA hand joints mapping: q r = q r 0  if  t = 0 ∫ t q ̇ r τ d τ  if  t ≠ 0 (31) with q ̇ r = J r w † q r J r s p r K s J h s † p h I Ω p h ̇ + K c p ̃ r − p r + J r w I − J r s K s J h s I Ω † J r s K s J h s I Ω u (32)

The camera system LeapMotion (Ultraleap inc., United States) is used to acquire active points of human hand. This camera system is composed of two monochromatic IR cameras and three LEDs that observe a hemispherical area and capture snapshots of up to 200 frames per second, and it returns information about positions of some fixed notable points of acquired hand as tips of fingers, the center of palm and knuckles. A ROS node acquires these from LeapMotion and publishes the values of active points, selecting only those useful for the mapping algorithm. To test the control law a virtual model of MUSHA hand is imported in V-rep (Coppelia Robotics GmbH, CH), and here its kinematic and dynamic behavior is simulated.

During surgery, the grasping and manipulation of tissue may generate both internal and external forces. The former impresses deformation and allows to avoid the sliding of the manipulated tissue, the latter generates twists. In accordance with the grasping theory (Prattichizzo and Trinkle, 2016), it is possible to evaluate from the applied force (i.e., λ) the components that generate the squeezing (i.e., ξ), and the ones that generate twisting (i.e., ω). This can be done by exploiting the grasping matrix G ∈ R 6 × 18 and the null space projector N G ∈ R 18 × 18 , following the relationship λ = − G † ⁡ ω + N G ξ (33) with G T = H G ̃ T . (34)

G ̃ T is the transposed grasping matrix that maps sources of twist from contact point reference to world reference. H is a selection matrix depending on the contact finger model. Considering the consistency of organic tissues, it is possible to approximate the grasp with the Soft Finger (SF) contact model. Interested reader is referred to Prattichizzo and Trinkle (2016) for further details. Therefore, admitting the applicability of the principle of kineto-static duality, it is possible to apply generalized forces exclusively in the directions of constrained motions. Directions of force application are shown in Figure 7 in case of contact with one finger. We assumed SF contact model for each point, thus the selection matrix H can be defined as H = H T 0 0 0 H I 0 0 0 H M ∈ R 18 × 18 (35) with H T = H I = H M = 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 . (36)

Consequently, the complete grasp matrix G ^T can be modeled as G T = H G ̃ T , being G ̃ T = R T ̄ T 0 0 0 R I ̄ T 0 0 0 R M ̄ T P T T P I T P M T , G ̃ T ∈ R 18 × 6 (37) with R ̄ j T = R j T 0 0 R j T   ∀ j ∈ T , I , M , (38) and P j T = I − S c j − p t s 0 I   ∀ j ∈ T , I , M (39) where R ^T ∈ R ^3×3 are the rotation matrices describing the orientation of each contact point reference frame Σ _j , j = T, I, M with respect to the world reference frame Σ _O , S(c _j − p _ts ) is the cross-product matrix of the vector c j p t s ̄ , and p _ts is center of mass of the grasped portion of tissue. Accounting the available grasping configurations, we can refine the general formula as follows: λ = N G ξ if method = P o − G † ω + N G ξ if method = P r − G † ω if method = R t (40)

At each time instant, the position of the centroid of the manipulated tissues and the positions of the contact points may be known, so the grasping matrix can be computed and a projector in its null space can be evaluated: N G = I − G † G . (41)

Thus, given measure of contact force, λ ̃ , the value of ξ and ω can be evaluated as: ξ = N G − 1 λ ̃ if method = P o N G † λ ̃ if method = P r 0 if method = R t (42) ω = − G † † λ ̃ if method ≠ P o 0 if method = P o (43)

The stability of the proposed control law was assessed by Lyapunov conditions. To this aim, we choose as positive definite Lyapunov functions the function of error: V = 1 2 ϵ T K c T ϵ . (44) Where ϵ is the distance between the actual (p _r ) and desired ( p r ̃ ) MUSHA fingertip. ϵ = p ̃ r − p r (45)

It is possible to evaluate its time derivative value as: ϵ ̇ = p r ̃ ̇ − J r w ̂ q r ̇ (46)

Consequently, the time derivative of the Lyapunov function is V ̇ = ϵ T K c T p r ̃ ̇ − J r w ̂ q r ̇ . (47) Where J r w ̂ is the Jacobian matrix of model of MUSHA hand evaluated in real current position. Substituting (28) in (47), the derivative Lyapunov function becomes V ̇ = ϵ T K c T p r ̃ ̇ − ϵ T K c T J r w ̂ J r w † M p h ̇ + J r w N M u + K c p ̃ r − p r . (48) Where M is the mapping action of control law M = J r s p r K s J h s † p h I Ω , (49) and N M is the projection action in null space of M used to overcome the duality of thumb pose in retractor method: N M = I − J r s p r K s J h s † p h I Ω † J r s p r K s J h s † p h I Ω . (50)

The desired velocity p ̃ ̇ r is chosen as the result of mapping M ( q r , p r , p h ) and projection in null space of mapping N M , p ̃ ̇ r = M p ̇ h + J r w N M u . (51)

Introducing (48) in (51), it is obtained that: V ̇ = ϵ T K c T p r ̃ ̇ − ϵ T K c T J r w ̂ J r w † p r ̃ ̇ + K c p ̃ r − p r . (52)

Reintroducing (52) the error evaluated in accordance with (45), we obtain V ̇ = ϵ T K c T p r ̃ ̇ − ϵ T K c T J r w ̂ J r w † p r ̃ ̇ + K c ϵ . (53)

If MUSHA is accurately calibrated the Jacobian matrix evaluated on the real position of fingers and on the estimated one are approximately equal, J r w ̂ ≈ J r w , and the time derivative Lyapunov function becomes V ̇ ≈ − ϵ T K c T K c ϵ . (54)

Since the Lyapunov function (47) is positive definite (V(t) > 0) and its time derivative is negative definite ( V ̇ ( t ) < 0 ) , the system is Lyapunov stable.

It is possible assert the BIBO stability, too. In fact, substituting definition of desired position (51) in control law (28), we obtain that q ̇ r = J r w † p ̃ ̇ r − K c ϵ (55)

It is possible to evaluate the tendency of error, using the definition of time derivative of error metric: ϵ ̇ = K c ϵ (56)

The error evolves as first order linear system, so the Lyapunov stability implies the BIBO stability.

Ten inexperienced participants (six males and four females, aged 24–59, all right-handed) took part in user studies. Each subject gave her/his consent to participate and was able to discontinue participation at any time during experiments. The experimental evaluation protocols followed the declaration of Helsinki, and there was no risk of harmful effects on subjects’ health. Participants performed 10 pick-and-place tasks of two cubes with 10 different initial scenarios pseudo-randomly sorted. The experimental set-up is shown in ¹ Figure 8. The dimension of each cube is 0.02 m × 0.02 m×0.02 m and the weight is equal to 20 g. To help the participant to complete the task, force feedback is provided by two vibromotors NFP-C1234L to inform them about internal and external forces applied on cubes. The two components of force are evaluated starting from the results of Section 3.1 using (42) and (43), respectively, and their norms in ℓ ² space, ‖ξ‖ and ‖ω‖, are used to calculate the motor PWMs as: P W M ξ = 1 k ‖ ξ ‖ Δ F 255 , (57) P W M ω = ‖ ω ‖ Δ F 255 . (58) where ΔF is the range of forces that can be acquired through MUSHA hand force sensors (ΔF = 4N). The dimension disparity between the internal force space ( ξ ∈ R 18 ) and the external one ( ω ∈ R 6 ) introduces an imbalance of haptic stimulus intensity which has been compensated by scaling the ℓ ² norm ‖ξ‖ with a k–factor equal to the ratio of the dimensions of the two spaces (k = 3).

At the end of each session, a NASA Task Load Index questionnaire (TLX) (Hart, 2006) was proposed to the participants, with the aim of assessing the perceived load in terms of Mental Demand (MD), Temporal Demand (TD), Physical Demand (PD), Performance (PE), Effort (EF), and Frustration (FR). This method assesses workload on a point-based scale. Each of the six questions has a scale of 21 levels, considering 1 as “very low” and 21 as “very high.” Questions are reported in what follows:

- Mental demand: “How mentally demanding was the task?”

All participants successfully completed the ten pick and place repetitions of the two cubes, so the completion times for each trial and the results of the questionnaires were considered as metrics for evaluation. Results are reported in Table 1 and Table 2. The mean time to complete the task is 16.49 s and no participant took more than 30s. The load of the task performed with the proposed control law is evaluated from answers of the NASA TLX questionnaire and its mean value is equal to 35%. According with the proposed task and experimental set-up we used the following weights: w _MD = 3, w _PD = 0, w _TD = 5, w _PE = 1, w _FR = 3, w _EF = 3. The physical demand was excluded from task load evaluation (PD = 0) because the physical load is due to the absence of arm support that is provided by the surgeon console in a real scenario.