Figure 3 (Left) shows a schematic cross section view of a typical PHRII. The human reference structure (the skeletal bone) is shown in the middle, completely encased by soft tissue both above and below it. The bone in this figure is the proximal phalanx of the index finger, and it's dorsal surface (above the bone in this figure) can be approximated to a straight line. Over the bone is a layer of soft tissue and skin represented by arrows showing its mechanical compliance in compression and sliding. Above the finger is the rigid plate of the robot attachment represented by a solid line, and its compliant padding in contact with the skin on the dorsal surface is similarly represented by arrows showing its mechanical compliance.

The PHRII region is defined by this contact surface described above, and it extends from the rigid reference geometry of the human to the rigid component of the robot. Also visible in Figure 3 (Left) is the antagonistic strap that is in contact with the skin on the palmar surface of the hand (below the bone in this figure). Identical to the considerations on the dorsal surface discussed previously, in this PHRII, we also have the two mechanically compliant regions of human soft tissue and device padding between the bone and the inelastic fibers of the strap which we consider as its reference geometry.

The Figure 3 (Left) also shows individual sources of movement as colored arrows. The strain distribution in the human soft tissue and the device padding are shown in purple and blue, respectively. Shear strain in the padding is also accounted for in our model, but not explicitly highlighted here. The human skin can also glide over the underlying layers of fascia, as shown in red. The reaction forces to this movement come from two main sources: stretching of the skin itself, which is dominated by the viscoelastic properties of the dermis, and the connective tissue linking the dermis to the underlying structures in the hypodermis and fascia. Lastly, we also have surface effects of slip and lift-off at the actual interface surface, as shown in green.

All components in the PHRII unit are modeled as viscoelastic elements. We use lumped parameter models to represent the robot and strap structures as linked point masses, and the human bone as the reference ground. With this approach, we represent the soft structures in the PHRII as viscoelastic elements connected in series as shown in Figure 3 (Right). This system can be viewed as being composed of three PHRII units, represented by the numbered circles. Each PHRII unit is connected to the human reference ground through series viscoelastic elements, and to each other either with compliant, or with inelastic constraints.

Figure 4 shows a bond graph representation of a single PHRII unit. The left half represents features in the tangential or x-direction while the right represents normal or the y-direction. Therefore, C_1y & C_2y represent the compressive elastic properties of the soft tissue and padding, while C_1x & C_2x represent the stretching/sliding of skin over layers of fascia and the padding shear strain. R_1x, R_1y, R_2x, and R_2y represent the viscous properties of each of these elements. The normal and the tangential strains are orthogonal and considered independently. However, they are still connected to the same inertia and reference ground. This individual unit is replicated at every PHRII and is used as the building block for the simulation environment. As seen in Figure 4 (Right), these units are coupled by both elastic and inelastic constraints, as well as external forces and this is the input to the primary lumped parameter model. The lumped parameter model is solved using a constrained ODE solved detailed in the next sub-section. The steady state and dynamic responses of the solver are the outputs which we use for further analysis.

In our example system, the dorsum of the hand has two such PHRII units connected by an inelastic constraint (simulating a rigid body) The palmar surface has a single one which is connected to the dorsal link via another viscoelastic element (simulating the strap). We will be discussing the construction of more complex systems with this building block as we go forward. Note that we have made the assumption of coupling the elastic elements in series separately from the viscous ones. Condition checks were included in the simulation of the system to ensure that this does not create unrealistic behavior in exception cases. These include negative normal forces at the interface on rapid withdrawal of loading.

We use experimentally obtained values for human tissue parameters where available in the literature. These include tensile viscoelastic properties of the skin (Silver et al., 2001; van Kuilenburg et al., 2012; Dabrowska et al., 2016), and compressive viscoelastic properties of the skin and the sub-cutaneous tissue (Wu et al., 2007). The thickness and size of all these regions across population samples are obtained from anthropometric data. The geometry and mechanical properties of the padding and straps can either be selected from a list of commonly used configurations or else optimized by the simulation to be returned as a design input.

Using the PHRII units defined above as our building block, we create larger systems for simulation. Figure 5 shows the system of the Maestro hand exoskeleton with our simulation model superposed over it. We first draw attention to the numbered point masses shown connected by inelastic constraints. Ternary links or sets of three connected masses (such as 1-2-5 and 6-8-7) move together essentially as a rigid body. The gray connections are elastic constraints to simulate compliance of the straps or in the exoskeleton linkage. Note that the human reference structure (shown as white linkages) will not be explicitly visualized in this system going forward, since it is considered as the reference ground. The proximal phalanx is considered as a moving reference ground based on transformations on the human MCP joint. Lastly, not all of the numbered points are in actual contact with the human (for example, points 5 & 8). The PHRII viscoelastic properties for these are set to zero while those for all others are set based on the human and padding viscoelastic properties discussed above.

These lumped parameter models are used in our dynamic simulations, allowing a system approximation that can be tuned to give reasonably accurate results while still being light enough to run in a real-time control loop. Since many of the experimentally determined properties of skin in literature have been approximated to piece-wise linear models (Silver, et al. Wu, et al.), our simulation also uses smoothed piece-wise linear models for tissue viscoelastic properties wherever possible based upon the desired complexity.

Dynamic simulations of the resulting multi-body system are formulated with ODE (Ordinary Differential Equation) solvers, and constraints are implemented with the option of simulating in both 2D and 3D using Lagrangian mechanics by Udwadia and Kalaba's method (Udwadia and Kalaba, 2002) as shown in Equation 1 (Here q's are the coordinates of each point, M is the overall inertia matrix, Q is the array of non-constraint forces for all points, and A & b are coefficients of the constraint equations after differentiating twice and grouping in the form of Equation 2). This method was chosen over penalty methods and simultaneous solutions for numerical stability and computational simplicity (Bayo and Ledesma, 1996; Witkin and Witkin, 2001).

To simulate dynamics through time, the ODE's can be solved both by implicit, as well as explicit methods. Standard ODE solvers are available for both kinds, with each having their respective advantages. Explicit adaptive step Range Kutta solvers such as Matlab's ODE45 are the most common type employed for differential equations. However, with numerically stiff systems, they tend to decrease step size drastically resulting in long computation times. Implicit solvers on the other hand are more resistant to this phenomenon and are capable of computing the dynamics of the system with far fewer temporal steps. However, care must be taken while using them as they are much more prone to numerical instability depending on the system parameters. The simulation environment makes use of both types of solvers to allow for explicit solvers for computational in general cases and for the use of implicit solvers in numerically stiff regions of the system.

The entire simulation was custom built in the Matlab environment (Mathworks Inc.) and provides both steady state and transient responses. However, due to the slow and almost quasi-static movement speed that is used in rehabilitation, assistance, and assessment, our results focused largely on the steady state response.

We configured HuRoSim to study the effects of applying an isometric extension load across the MCP joint. The input load to the simulation is in the form of an applied moment from the exoskeleton (to a maximum of 0.5 Nm), and the studied outputs are the final displacements and forces present in the system on reaching a steady state.

The human and exoskeleton interface material parameters were gathered from literature sources, most notably (Silver et al., 2001) which establishes the piece-wise linear approximation of stiffness properties of human skin caused by the behavior of different fibers (elastin and collagen) within the dermis of the skin. In the tangential direction, the values from this source start with an initial elastic stiffness of 0.1 MPa upto a strain of 0.4, followed by a sharp increase trending up to 18.8 MPa. Wu et al. (2007) similarly establishes values of compressive stress-strain curves for both skin and sub-cutaneous tissue that were used in the simulation. It should be noted that detailed in-vivo measurements of the tissue properties of the surfaces of the human hand have not been widely studied, leading to a lack of available literature to obtain measurements from. This scarcity contributed to the motivation for indentation experiments in section 2.3.3 that gather in-vivo stiffness data of the hand dorsum for use by HuRoSim.

The device parameters were measured directly on the maestro hand exoskeleton, having a dorsal rigid surface of 50 mm length for the metacarpal region and 30 mm for the proximal phalanx of the index finger, with strap widths of 10 mm.

From the applied moment loads to the isometric finger and the tissue properties chosen, we are able to simulate the dynamic loads that develop. These load patterns generate kinematic behavior that we then tested with an experimental set-up with the Maestro hand exoskeleton and the hand.

We test the predicted kinematic behavior from the simulation with an experimental test bed. We instrument the maestro hand exoskeleton Figure 6 to enable easier measurement of the effects of applying an isometric load across the MCP joint. The index finger was locked about the MCP, PIP, and DIP joints, and an external moment identical to that applied in the simulation was developed over the hand.

Only the index finger of Maestro was attached to the human test subject, and loads were applied via Maestro's MCP joint. The subject's arm is immobilized at the forearm and wrist, in addition to the palmar surfaces of the hand, fingers, and thumb. Phase space X2E motion capture markers were attached to multiple points on key linkages of the maestro exoskeleton. Markers were also placed on the metacarpal, MCP and PIP joints of the index finger and the base of the thumb to track the position of the palm and the index finger. This enables a fair comparison to the simulated system by confirming that the movement of the human hand is minimal.

The required moment load was applied using the Maestro exoskeleton's torque sensing series elastic actuators. To achieve a known level of strap tension, pre-tensioning was carried out using a calibrated Omega DFG55 load cell, rated for 25N. Strap pre-tensioning was performed prior to placing the hand in the brace of the experimental setup by affixing the Maestro grounding link to a rigid fixture through a dovetail joint and aligning the hand within it. The straps were then pulled to the target tension with the load cell attached in series. To ensure minimal effects of friction between the strap and the palm during this step, both ends of all straps were pre-tensioned individually in a consistent order of: distal-lateral, proximal-lateral, distal-medial, proximal-medial. A maximum displayed error of 0.1N was allowed in the experiment protocol during the fastening process.

To design an optimal PHRII for the hand dorsum, we simplify the complex interaction at the interface as two plates of length L and uniform width sandwiched between the rigid reference plate of the Maestro robot, and the rigid human bone. The robot reference plate is held to the dorsum by a bias force mimicking a strap (F_b), applied normally and at L/2 as seen in Figure 8 (Left).

We then load the PHRII with an external force mimicking the reaction forces from the Maestro actuators. Using the principle of transmissibility, we express the applied external force to the attachment plate as a combination of an equivalent force and moment applied at the center of the plate, placed coaxially with the bias force. The resulting reaction pressure (P_reaction) distribution balances the net force (F_b plus the normal component of F_L) and the external moment (M_L) [Figure 8 (Right)].

Our objective is to balance the applied force and moment while minimizing peak reaction pressure (C_i) along the contact surface between the two plates.

Under the counter-clockwise external moment, M_L, applied on the system, minimizing the cost function, C_i gives us a pressure distribution with two regions of optimized uniform reaction pressure (P_opt) below the plate (Figure 9).

Since we assume that the sum of pressure is equal to the applied load, minimizing the peak pressure would distribute it over a larger area giving us a region of uniform pressure. However, uniform pressure across the entire plate cannot balance the counter-clockwise applied moment, which explains the existence of two sections of different uniform pressure. The direction of this reaction pressure distribution (seen in Figure 9) depends on the direction of the external loading and reverses itself when the external loading is in the clockwise direction.

The boundaries of the regions of uniform pressure shift depending on the ratio of external applied moment (M_L) to the external applied force (F_L). The position of “x,” the center of the highest pressure region in the optimal distribution (P_opt), and “P_peak,” the magnitude of the highest pressure between the two uniform distributions, are computed in Equations (2) and (3).

The optimal value of bias force, F_b of the attachment against the dorsum, to minimize C_i for the given force and moment loading configuration, is calculated to be the lowest value that gives us a non-negative pressure region (4):

By analyzing the possible solutions to achieve this optimal pressure, (P_opt), a hyperbolic distribution with high stiffness in the center of the plate which tapers off toward the edge in the direction of the applied moment, M_L, is a solution to the equation, though its exact parameters are yet to be computed. We approximate this hyperbolic representation to a linear stiffness profile for our experiment, decreasing from the highest stiffness at the middle of the attachment (k_mid), and symmetrically tapering off to a minimum at each edge (k_edge). The symmetry allows the resulting profile to hold true for external moment loads in either direction. To obtain k_pHRI, we need to obtain k_dorsum in addition to knowing k_eff.

We used the HuRoSim environment to compute the effective spatial stiffness gradient that minimizes peak pressure. The environment was also used to characterize the relationship between the bias force, the gradient of stiffness, relative displacement between the robot and the human, and the peak pressure over the hand dorsum.

For this experiment, the dorsum surface was discretized into 15 total points, with the PHRII interacting with the underlying human metacarpal through the stiffnesses k_pHRI and k_dorsum in series at each point. Piece-wise linear values for both stiffness were used, with k_dorsum taken from our experimental results. The system was simulated for varying applied force (F_L) and moment loads (M_L) with varying k_pHRI profiles to examine the resulting pressure distribution. The relative displacement of the PHRII with respect to the underlying bone, due to F_L, was also captured for each stiffness profile, and these results are presented in the next section.

Tissue stiffness on the human body is quantified by applying force over a range of displacement.

To measure k_dorsum, we designed an indentation system comprised of a Phantom Premium 1.5 high force haptic renderer, which has a high positional accuracy of the end effector (7*10⁻⁶m). This was used along with an ATI Nano 17 force torque transducer (having a high force torque sensing accuracy of 0.001N) attached at the end of the linkage as an indenter to probe the hand dorsum (Figure 10). We selected five points along a line between the metacarpophalangeal joint and the radial styloid process along the 2^nd metacarpal bone, the 3^rd metacarpal bone, and along a line between the two metacarpal bones, in the inter-metacarpal region. This region was selected to correspond to the area of the attachment plate on the Maestro exoskeleton. Figure 10 (Middle) shows the regions selected for indentation.

The wrist and arm were supported in braces, and a consistent grasp object was used across trials to minimize its influence on stiffness distributions due to the changes in bone locations and muscle recruitment strategy, as seen in Figure 10 (Left). The wrist and arm supports were positioned to level the hand dorsum in the transverse anatomical plane. The phantom probe is then manually led once to each point marked on the dorsum for indentation. The probe uses these points as input into an interpolator to compute a spatial trajectory to follow. The phantom is driven in an open-loop position controlled configuration.

Quantifying the stiffness of the indentation system in the direction normal to the hand dorsum demonstrated the need to account for this value in estimating the stiffness of the hand dorsum. The stiffness of the indentation system (k_indenter) was found to be 2.67 N/mm along the workspace. We account for this stiffness in the measurement of the hand dorsum stiffness (k_dorsum) by modeling the interaction between the indentation system and the hand dorsum as two springs in series. The measured stiffness of the hand dorsum, k_measured, (Figure 10 Right) is used along with k_indenter to calculate the k_dorsum Equation (5).

We made five sequential repeated measures to estimate the variance in the measured stiffness at each point. The observed variance was attributed to movement in the hand. The hand dorsum stiffness data is collected from one pilot subject only, and with a probe having a square base with 4mm edges and 1.5mm filets on each edge to minimize discomfort during indentation. The indenter profile and size were chosen iteratively based on the relative distribution of the hard and soft tissue structures in the hand. Increasing the resolution of the grid beyond the current levels introduced errors due to partial overlap of soft and hard tissue at the points of measurement.

In our model, k_eff between the human bone and the Maestro robot's reference plate is comprised of k_dorsum and k_pHRI in series. Therefore, once we have numerically computed k_eff, and measured k_dorsum through the indentation experiment, the required k_pHRI can be calculated at every point on the attachment surface Equation (6). This gives us a PHRII stiffness profile that should generate the minimum peak pressure, or the optimal pressure profile (P_opt) on the hand dorsum for the given F_L and M_L.

The effect of varying the applied moment loading and the strap pretension were studied in both the HuRoSim environment, as well as in the physical experiment. Figure 11 (Left) shows the angular displacement of the exoskeleton dorsum attachment on application of moment load by the exoskeleton about the MCP joint. HuRoSim predicted an angular displacement that increased approximately proportional to the applied moment load, and that was reduced on increased strap pretension. Physical experiments under the same conditions confirmed both trends as seen in Figure 11 (Right). The experimental values of angular displacement under low strap pretension matched the predicted values quite closely. However, under higher pretension, where the non-linear behavior of human soft tissue becomes more pronounced, there is expected divergence between the predicted and experimental values.

In addition, when the exoskeleton applies an extension torque about the MCP joint, HuRoSim also predicted an increase in pressure at the location on the PHRII on the hand dorsum closer to the joint while the location further away from the MCP joint experiences a reduction in pressure. The inverse is true when the actuator applies a flexion torque about the MCP joint. HuRoSim indicated that the magnitude of these increases with respect to the applied moment was found to be dependent on the length of the exoskeleton attachment and on the geometry of the PHRII at the contact locations on the hand dorsum.

As a consequence of the distribution of pressures described above, HuRoSim predicted that the end of the PHRII closest to the MCP joint would lift off and lose contact with the skin surface on the application of extension torques by the exoskeleton. When we replicated the simulation empirically with the Maestro hand exoskeleton, we observed the lift off phenomenon predicted by HuRoSim (Figure 12).

The effect of varying strap stiffness was studied using the simulation environment to identify the optimal strap stiffness to minimize both output metrics (relative displacement and forces across the interface). HuRoSim's initial results indicated that the best performance on both metrics occurred with effectively inelastic straps, or straps with extremely high stiffness.

This result was unexpected since pilot user feedback with the maestro exoskeleton reported that mild elasticity was the most comfortable. Exploring different variables in the system, it was observed that in the actual human hand, co-contraction of muscles in response to the applied load leads to an increase in the thickness of the muscle body as it is flexed, and resulting changes in the cross sectional thickness of the hand. These also vary depending on the musculature of the individual subject and the magnitude of voluntary contraction. When the simulation was tested with approximate values of the cross section variation, a slightly reduced strap stiffness was found to produce lower pressure distributions with minimal trade-off in the relative displacement values.

The values of this optimal reduced stiffness, however, are dependent on the magnitude of cross section increase due to co-contraction. Future work can perform the simulation using magnitudes of cross sectional variation sourced from literature or experiment to inform design values of this optimal stiffness.