used in many areas, such as medicine, rehabilitation, industrial environments, and augmented reality platforms [1], [2], [3], [4], [5]. Beyond supplementing human motor capability (such as when using exoskeletons), allowing kinesthetic haptic feedback in complex motor tasks can create an augmented perception that can improve outcomes compared to systems without haptics [6]. However, there are some challenges in designing controllers for haptics-enabled robotic systems and, in general, when humans and robots exchange physical energy through continuous kinesthetic contact. The trade-off between stability and transparency is known as a central control challenge for physical human-robot interaction (pHRI). This means that in order to have an ideally perfectly transparent system (which would have an ideally high resolution for rendering very low to very high contact impedances), the closed-loop system would be pushed toward the edge of instability. A perfectly transparent system can become unstable in the presence of the smallest delay and possibly noise and uncertainties [7], [8]. Thus, most pHRI controllers adopt a conservative control solution that sacrifices the quality of interaction to guarantee the stability of the humanrobot interaction. Stability must be met due to the close contact of the powered robot mechanics and human biomechanics. In some cases (such as exoskeletons or rehabilitation robots), humans are rigidly attached to the body of the robot, and any out-of-control instability could severely damage the engaged joints of the involved human.

There are several established methods for guaranteeing stability, such as the widely-used Time-Domain Passivity Approach [9], [10], which works by injecting damping into the system when non-passive energy is detected. However, the classic derivation of this method can be too conservative and hinders haptic transparency. As a result, a wide range of extensions to this method have been proposed in recent years that attempt to reduce the conservatism and thus increase haptic transparency [8], [11], [12], [13], [14]. However, many of these approaches overlook an important aspect of human-robot interaction -the ability of the human limb to absorb and dissipate some energy from the system through the reactive component of the biomechanics.

We have recently shown that this capability [15], [16], [17], [18], quantified as the Biomechanical Passivity Index (BioPI), can be used as a passivity margin in the synthesis of passivitybased stabilizers to allow some non-passive energy to pass through to the user while still guaranteeing stability. This is necessary for some pHRI scenarios, such as assistive robotic physical therapy or assistive exoskeletons. In these situations, non-passive energy should be delivered to the user to assist with their desired motion [15], [19], [20]. In traditional TDPA, this energy would be damped out, and the user would not be assisted [15], [16]. However, utilizing the BioPI as a stability margin allows some non-passive energy to pass to the user, allowing them to be assisted without risking instability. The example control designs based on the assumption of having access to a measure of BioPI can be found in our recent publications [15], [16]. However, to maximize the performance of BioPI-based passivity stabilizers, a less conservative and more accurate estimation of BioPI would be needed.

In order to assess the biomechanics of the limb, an offline identification process must take place, where the reactive dynamics of the user are read from the robot's torque and velocity sensors and used to create a model. This approach is common across different methods of assessing biomechanics, including linear methods [21], [22], [23], such as impedance and massspring-damper modeling, and non-linear methods [16], [24], such as passivity. The advantage of using passivity is that the non-linear nature of the limb biomechanics can be represented. In this work, we go beyond the offline identification process to create a mapping between the muscle activity and BioPI that can be used in real pHRI tasks, even when the torque and velocity readings are superimposed with the actions of the user. Even a conservative estimate of BioPI can significantly enhance the transparency of a robotic system [15], [16], [25].

It has previously been shown that the BioPI varies between users and conditions, such as the direction of interaction and the grip strength of the user [16], [17], [18], [26]. Furthermore, we have recently shown that there is a correlation between the magnitude of surface electromyography (sEMG) signals from sensors placed on the forearm of the user and their BioPI [17], which motivates the generation of models based on this biosignal modality. However, the relationship between BioPI and muscle activity is more complex than can be captured by only utilizing the magnitude of sEMG signal from individual sensors. This is because the intricate co-contractions and co-modulations of muscles in the upper-limb affect the viscoelasticity of the biomechanics and, thus, the energetic behavior. So, there is a need to design a computational model that can look into the activation of various muscles holistically and in one frame to capture the effect of co-modulation and decode the corresponding change on the passivity of the reactive dynamics of the limb. To better capture this complexity, in this work, we propose, for the first time, to utilize the emerging concept of a 'functional muscle network' to create a model for real-time estimation of BioPI.

Muscle networks use concepts from network theory to encode information about the connectivity between muscles in the body. Each sEMG sensor represents a node in a graph, and the connectivity (which can encode the coupling) between each sensor pair determines the weights of the graph's edges. In this case, magnitude-squared coherence was used as the connectivity metric. This method has previously been used to enhance understating of human biomechanics, such as identifying biomarkers for disease and analysis of muscle coordination during tasks [27], [28], [29].

In this work, we construct muscle networks using signals from an array of sixteen sEMG sensors placed on the forearm and use the edge weights to build regression models to predict the BioPI of the user. Additionally, it is desirable to limit the number of sensors and thus generate smaller graphs needed for this estimation, with the goal of reducing the computational cost and setup time. As a result, we analyze different sizes of subnetworks to determine how the performance of the BioPI prediction changes with the network size and if an optimal network can be proposed to secure a high-quality prediction with a minimal graph size. It should be highlighted that the generalizability of any human-centric model for pHRI is always an area of interest. Thus, we conduct a generalizability analysis in this study to evaluate differences in performance when using an optimal subnetwork specific to each user, compared to a generalized subnetwork that is optimized for the best group performance. Thus, in summary, the specific goals of this work are:

r Determine if the functional muscle connectivity informa- tion contained in muscle networks can be used in a regression model to accurately predict BioPI.

r Investigate the minimum subnetwork size to predict the BioPI without compromising performance.

r Evaluate if there is a performance cost to selecting the nodes in the subnetwork at the group-level compared to the subject-level.

r Explore patterns in the selected nodes of the optimal sub- networks.

In this work, we consider the interconnected system composed of a human (physically interacting with the robot) and the robotic system. It has been shown that such a system is guaranteed to be L2 stable if it meets the requirements for an Output Strictly Passive (OSP) system, as defined in non-linear control theory [18], [30]. The definition of OSP is given by:

In (1), U (t) and Y (t) are the input and output vectors of the system at time t, respectively. E(t 1 ) is the initial energy in the system and T 1 and T 2 are the start and end times of the interaction. ξ represents the BioPI of the system. For ξ ≥ 0 (necessary for an OSP model), this is also known as the Excess of Passivity. OSP systems are L2 stable, with an L2 gain of 1/ξ. In contrast, if ξ < 0, the system is Output Non-Passive, and ξ represents the Shortage of Passivity.

Utilizing this definition of OSP, the BioPI of the human biomechanics can be computed via system identification. For this measurement, the user holds onto the robot while it perturbates the limb. By recording the interactional torques and angular velocities, the true BioPI ( ξ) of the user can be calculated as follows:

where τ is the applied torque (input to the system) and ω is the angular velocity (output of the system). More details on this approach can be found in [15], [16], [17].

As previously noted, knowledge of the real-time value of the BioPI during a robotic task can be used to improve the performance of a robotic controller and allow for enhanced haptic feedback [15], [16], [25]. However, this system identification method of computing BioPI is only feasible when the user is not actively moving the robot. In active pHRI tasks, the motions and forces of the user would be superimposed on the information needed to identify the BioPI. Thus, there is a need for an additional source of information. In this work, for the first time, we proposed the use of a muscle network model (based on multichannel sEMG recordings) to estimate the BioPI. Fig. 1 shows a proposed use case for this information. The dashed lines represent the new contributions from this work compared to traditional TDPA, where the BioPI estimate is assumed to be zero.

This study includes ten human subjects (five male, five female, mean age 27.7 ± 6.8 years). The Institutional Review Board of New York University approved this study. All subjects signed a written consent prior to participating in the study. Subjects denied any history of musculoskeletal injury.

An overview of the setup and procedure is shown in Fig. 2. During the experiment, each subject gripped onto the handle of a High-Definition Haptic Device (Quansar, Markham ON, Canada) robot while it perturbed their wrist by rotating in one of four directions: abduction, adduction, pronation, and supination. All subjects were right-handed and held the robot handle in their right hand. The perturbation signal was a mixed-harmonic sine wave with frequencies ranging from one to five Hertz (covering the typical range of voluntary human interaction frequencies [31], [32]). The perturbation angle was centered at fifteen degrees from vertical and varied from zero to thirty degrees. Each trial was one hundred seconds in length. This length was chosen as subjects had difficulty maintaining a consistent, stiff grip for longer periods of time. Of the ten subjects, four were experienced users of the robot, and the remaining six completed two to four practice trials to acclimatize to the robot.

Subjects wore sixteen Bipolar Delsys Trigno sEMG sensors (Delsys, Natick, MA, USA) on their forearm throughout the experiment. The sensors were placed at four distances from the wrist, along the Brachioradialis (Sensors 1-4), Flexor Carpi Radialis (5-8), Extensor Carpi Ulnaris (9)(10)(11)(12), and Extensor Digitorum (13-16) muscles [33]. Bipolar sEMG sensors scanning the arm were chosen over high-density sEMG scanning a smaller area because the focus of this work is on the interconnectivity between muscles. Prior to starting the experiment, the Maximum Voluntary Contraction (MVC) of the subject was recorded for two sensors (Sensors 8 and 16). During each task, subjects were shown live visual myofeedback of their muscle activity from these two sensors. Two grip conditions were tested in this experiment: a stiff grip and a relaxed grip. For the stiff trials, subjects were instructed to maintain thirty percent of their MVC, and for the relaxed trials, subjects were instructed to maintain five percent of their MVC while keeping all their fingers and their palm in contact with the robot handle. These grip conditions were chosen to investigate the change in BioPI under different grip strengths. Thirty percent MVC was found to be the highest value that subjects could maintain consistently over time, while five percent MVC was the minimum to maintain contact with the robot handle. In all trials, the goal for the subject was to maintain a steady grip on the robot handle without exerting additional force on the robot by moving it (other than the natural response of their biomechanics).

For each perturbation direction and grip condition, two trials were completed, and these were randomly split into either the Test or Train data sets. The order of the trials was randomized, and a fifteen-minute break was taken after each stiff trial to prevent an accumulation of fatigue from affecting future trials.

The posture and positioning of the subjects were instructed to minimize variation between trials. Participants stood with the upper arm against their torso and elbow bent at a right angle, with their forearm extending directly ahead of them. The robot was placed on a height-adjustable table to allow all subjects to achieve this position. The subjects gripped the handle naturally, with the thumb on one side and fingers wrapping around the other side, to mimic the conditions in a pHRI task. Subjects were asked to remain still throughout the trial, and the position of the robot handle was held constant.

The data was processed in MATLAB (Mathworks, Natick, MA, USA). The BioPI data was epoched into one-second segments and averaged over each segment. The first five seconds of each trial were discarded due to transient artifacts in the computation when beginning the trials.

For each sensor and trial, the sEMG signals were filtered using a 4th Order Butterworth bandpass filter between 20 Hz and 500 Hz. Subsequently, 4th Order Butterworth bandstop filters of width 4 Hz were applied at integer multiples of 60 Hz to remove powerline noise. The sEMG signals were divided into one-second segments to match the BioPI data. The first five seconds of data were discarded for each trial. For each segment, the magnitude-squared coherence (MSC) between each pair of sensors was computed using MATLAB's mscohere function, with an epoch length of 64 points (approximately 32 milliseconds) and a 50% overlap. MSC is a metric of similarity between two signals in the frequency domain, which, in this context, is used to quantify the connectivity between two muscles. For two signals, X and Y, it is computed as:

where P xy is the cross-spectral density between signals X and Y, and P xx and P xx are the power spectral density of signals X and Y, respectively [29]. The mean MSC was taken over the range of 20 to 50 Hz for each sensor pair.

For each subject and trial, muscle networks were created at each time point using the sensors as the nodes and the MSC between each sensor pair as the weights of the edges in the graph. That is, the input to each muscle network was the sixteen sEMG signals, and the output was a graph of sixteen nodes with a weighted edge between each pair of nodes corresponding to the connectivity of the muscles at those sEMG locations. These edge weights were then used to predict the BioPI throughout the trials. We have recently shown that there are large differences in the BioPI in the abduction and adduction (A/A) directions compared to those in the pronation and supination (P/S) directions [17]. As a result, the analysis of the BioPI was conducted separately for A/A versus P/S.

The optimal subnetworks of each size, from two to sixteen nodes per network, were investigated. To find the bestperforming subnetwork with n vertices, every combination of n out of sixteen sensors was considered for each subject.

Using only the data from the Train data set, the edges of each resulting subnetwork were fed as the input into a regression model to predict the BioPI. Though more complex machine learning algorithms could also be employed here, this is outside of the scope of this work. The regression model was trained to find the best fit between the inputted edges and the true BioPI. This is done by minimizing the sum of the squared error between the estimated BioPI and the true BioPI. For each of these subnetwork models, the R 2 correlation between the predicted and true BioPI values was computed. The subnetwork that resulted in the highest R 2 value was selected as the optimal subnetwork for that network size and subject. It should be noted that this analysis is at the subject level; thus, the optimal subnetwork for subject i may have a different combination of sensors when compared to the optimal subnetwork for subject j. The resultant trained models were then evaluated using the Test data set for each user. The R 2 values between the predicted and true BioPI for each subject were used to assess the models and compare the impact of network size.

As the next phase of the study, instead of conducting the investigation at the subject level, we select the generalizable network that produces the highest mean R 2 value across all subjects. In other words, this subnetwork's composition performs the best among all subnetworks when the goal is to select the same network composition (sensor location and number) for all subjects.

A statistical analysis was performed to compare the individual and generalized subnetwork results. The Kolmogorov-Smirnov test rejected the null hypothesis that the data was normally distributed. Following this test, Wilcoxon signed-rank tests were completed to assess the performance differences between the individual and general subnetworks at each network size with a significance level of p < 0.05.

Following the steps in Section II-D, BioPI prediction models were produced for each subject and network size. The output of the models evaluated on the Test data sets are shown for the subjects with the median R 2 values between the Predicted and True BioPI in both direction sets in Fig. 3 (Subject 4 for A/A and Subject 1 for P/S). For both the directions shown, we see that the median participant has a high R 2 value for all the network sizes shown, with a modest improvement as the network size increases. For the full sixteen-sensor network, the median R 2 is 0.88 for A/A and 0.90 for P/S. This indicates that there is a high correlation between the muscle network edges and the BioPI.

In all the plots in Fig. 3, there is a clear separation between the True BioPI for the relaxed condition (low BioPI) and the stiff condition (high BioPI). In some cases, there is also separation in the True BioPI based on the direction of interaction. This can be seen in the A/A scatter plots, where there are two clusters of points in the high BioPI region, representing the Abduction and Adduction trials individually. In general, there is little variation in the True BioPI throughout a given trial. However, the Predicted BioPI can have greater variation, which results in the 'columns' of data points visible in the plots.

To analyze the effect of the subnetwork size on the model performance, consider Fig. 4. The purple boxes represent the R 2 values for each subject using the subnetworks individually optimized for each subject. The top plot shows the results for the A/A case, and the bottom plot shows the results for the P/S case. In both plots, the performance tends to improve as the number of sensors in the network increases. These performance increases are larger for the smaller network models (e.g., moving from a two-sensor network to a three-sensor network) and then level out for the larger network models. This indicates that the information gained from adding another sensor to a smaller network is larger than that from adding another sensor to a larger network. This is logical since a larger network may already contain much of the information included in the new node. We see in the boxplots that the performance levels out at approximately eight sensors in the A/A case and six sensors in the P/S case. Reducing the number of sensors helps with the practicality of this approach for real-world pHRI tasks by lowering the time and complexity of the setup.

Another relevant consideration here is that as the number of sensors in the network increases, so does the number of regressors in the model. With more regressors, a better fit can be obtained, but there is a risk of overfitting. Since the model is built using the Train data set and then applied to the Test data set, any overfitting is accounted for in this analysis. This explains why, in the A/A case, the median R 2 actually decreases from the thirteen-sensor model to the sixteen-sensor model. This result can also be seen in the P/S model, where the ten-sensor model outperforms the sixteen-sensor model.

We also compared the results of the individual subnetwork models to those of the general subnetwork models, as shown in Fig. 4. The individualized results are shown in purple, and the generalized results are shown in orange.

For both the A/A and P/S cases, there is a large performance drop between the individual and general models for the smaller subnetworks. However, for the larger subnetworks, the performance is more similar, and, in many cases, there is no statistically significant difference between the two models for larger subnetwork sizes. For the A/A case, performance evens out for subnetworks with at least seven sensors, and for the P/S case, performance evens for subnetworks with at least ten sensors. One explanation for this phenomenon is that there is more variation between subnetworks at smaller network sizes. For example, for the two-sensor subnetworks, there could be very little overlap between the optimal sensor selections across subjects. However, in larger networks, where more sensors must be selected, there is necessarily more overlap in the sensors selected by the two models. Hence, the general model is closer to the individual model. In the sixteen-sensor case, there is only one possible subject (selecting all sensors), so the results are identical. Thus, while there is a performance cost to using a general network, particularly for cases with low numbers of sensors in the network, that cost is minimal (and, in many cases, not statistically significant) for subnetworks containing larger numbers of sensors. This insight is useful for practical applications in pHRI. Using a generalized network reduces the need for recalibration, reducing the setup time and computational cost for a new participant to use the system.

Figs. 5 and 6 show visualizations of optimal subnetworks with four, six, and eight sensors. For Fig. 5, the results are shown from the Abduction Test trial at the 10-second time step. Similarly, in Fig. 6, the results are shown from the Pronation Test trial at the 10-second time step. For each figure, the results are shown for the subject with the median performance (as described in Section III-A). For the A/A case, this was Subject 4; for the P/S case, this was Subject 1. The thickness of each line indicates the MSC between the given sensor pair, with thicker lines denoting higher MSC; thus, higher coherence between the node activities (the sEMG signals).

Considering the four-sensor network in the A/A case, we can see that the MSC between the selected sensors (for this subject, Sensors 4, 8, 15, and 16) increases noticeably for the stiff grip compared to the relaxed grip. This trend is seen in all the connections for the four-sensor network and is particularly prominent in the connection between Sensors 4 and 8.

Next, considering the six-sensor subnetwork, it is notable that all of the sensors that were selected in the optimal four-sensor subnetwork are selected again in this subnetwork, with the addition of Sensors 9 and 11. This pattern of larger optimal networks containing the previously found smaller optimal networks was consistent across all subjects and in both the A/A and P/S cases. Comparing the edges of the six-sensor network for the different grip conditions, there is a similar trend here as in the four-sensor network. For example, for the edge between Sensors 4 and 11, MSC tends to increase for the stiff grip compared to the relaxed grip. However, some of the sensor pairs do not show much difference between the relaxed and stiff grip, such as the edge between Sensors 9 and 11, which has similar weight in both cases. This indicates that while some useful information has been added by increasing the size of the network, some of the other edges that were added may be less useful in distinguishing between grip conditions (and, thus, the BioPI).

We see a similar phenomenon when comparing the six-sensor subnetwork to the eight-sensor subnetwork. While some of the new edges show a large difference between the relaxed and stiff conditions (such as the edge between Sensors 5 and 16), other edges are very similar between the two conditions (such as the edge between Sensors 5 and 9). Thus, as the subnetworks get larger, they encompass more information about the biomechanics, which can be used to predict the BioPI. However, these larger networks also incorporate edges that may not contain information relevant to the BioPI prediction, which could cause the regression to overfit. These observations track with the results discussed in Section III-B. The marginal gain in model accuracy for adding a new node decreases in larger subnetworks, as most of the information is already captured in the network.

There is a similar trend in Fig. 6, where all the edges in the four-sensor network respond to the change in grip condition (with a particularly strong response for the edge connecting Sensors 3 and 15) while the larger networks have more diversity in response to the change in condition. Additionally, in the six-sensor case, we can see an example of the edge weight decreasing between the Relaxed and Stiff conditions in the edge connecting Sensors 4 and 8. Although it may seem unintuitive that the coherence would be lower while the muscle activation is increased, this situation highlights the lack of synchronicity between the two signals. This underscores the importance of using coherence and the holistic concept of muscle network rather than just looking into the absolute activation of each node in predicting the BioPI.

In comparing the results in Figs. 5 and 6, one can also note the differences in the networks between the subjects and directions. While the edge between Sensors 4 and 8 shows a large increase between the relaxed and stiff conditions in Subject 4 for the A/A direction, there is instead a decrease in the connectivity between those sensors for Subject 1 in the P/S case. This example highlights the usefulness of a subject-specific and direction-specific model for the BioPI prediction since the biomechanical response can vary across users and perturbation conditions.

Zooming out, we can also consider which sensors were selected most often across all subjects for subnetworks of different sizes. Fig. 7 shows the selection frequency for each sensor for the optimal subnetworks with four, six, and eight sensors. The results are shown separately for the A/A and P/S cases. It can be mentioned that, in general, there are variations between subjects regarding which sensors are selected. However, there are some important trends to be highlighted, which can indicate that some sensors are generally more useful than others. First, it should be Fig. 7. Each bar represents the frequency of selection of the sensors in the optimal subnetworks across subjects. The x-axis indicates the sensor number, and the y-axis indicates the number of subjects for which the given sensor was included in their optimal subnetwork, divided by the total number of subjects. The results are shown for networks of size four (top), six (middle), and eight (bottom). The A/A results are shown in orange, and the P/S results are shown in purple.

noted that for the four-sensor subnetworks in the A/A case, five of the sensors are selected for only one subject, and in the P/S case, four of the sensors are selected for at most one subject. For both directions, the most frequently selected sensor in the four-sensor networks was Sensor 8, which was selected for six of the ten subjects. Considering the six-sensor networks, we see that over half of the subjects selected Sensors 3, 8, and 15 in the A/A case. Similarly, over half the subjects selected Sensors 1, 8, and 10 in the P/S case. In the eight-sensor networks, at least seventy percent of the subjects selected Sensor 3, 8, 9, and 13 in the A/A case, and at least eighty percent selected Sensors 1, 5, and 8 in the P/S case. The presence of these highly selected sensors helps explains the high performance of the larger generalized subnetworks, as shown in Section III-C. Since many of the same sensors are being selected for the majority of subjects, the generalized subnetwork is similar to the individually selected subnetworks.

Since the models are built using connectivity between sensors instead of the sensors themselves, it is also interesting to consider which pairs of sensors are selected together most often. In Fig. 8, the frequency that sensor pairs are selected across subjects both the A/A and P/S cases. As previously, the results are shown for the optimal networks with four, six, and eight nodes. These heatmaps highlight the variations in the frequency of selections across sensor pairs.

Considering the four-sensor subnetwork results in both directions, we see that most sensor pairs (edges) are never selected together. There are some sensor pairs that appear in the optimal network for multiple subjects, but never more than three out of ten subjects. However, in the six-sensor networks, some sensor pairs are selected more frequently. In the A/A case, half of the subjects have edges (3,15) and (4, 15) selected. In the P/S case, edge (1,8) is selected in half of the subjects. Since only an eighth of all edges can be selected in a six-sensor network, it is notable that some edges are found in the optimal networks for so many subjects. Considering the eight-sensor networks, sixty percent of the optimal subnetworks in the A/A direction include edges (8,13) and (8,9). In the P/S direction, there is even more coordination, with eighty percent of the optimal subnetworks including edges (1,8) and (5,8) and seventy percent including edge (1,5).

Again, these commonalities between the optimal subnetworks, especially for the larger subnetworks, explain the high performance of the generalized networks. The overlap between the subjects means that the nodes selected for the generalized subnetwork can closely match the individually selected subnetworks.

IV. DISCUSSION AND CONCLUSION Knowledge regarding the energy absorption capacity of the human user's limb during interaction with a robotic system can reduce the conservatism necessary to maintain stability and thus allow for improved control of pHRI systems and haptics rendering. Previous work in this area has relied on offline system identification procedures, which cannot adapt to changes in the user's biomechanics. However, in this work, we evaluated the interplay between the functional muscle networks and the energetic behaviors of the upper-limb at two levels of co-contraction and four directions of perturbation to create an accurate model of the BioPI that can update in real-time during dynamic pHRI tasks. To reduce the setup time and computation complexity, we investigated the performance cost of considering only a subset of sensors (subnetworks) in the BioPI model, resulting in smaller graphs explaining distributed co-modulation between the nodes/muscles. Our investigation showed that even a subnetwork of six to eight sensors can perform similarly to the full sixteen-sensor network. Additionally, we considered the performance difference between a subnetwork optimized for each individual subject compared to a generalized subnetwork. For larger networks, the generalized and individual models had similar performance, which emphasizes the possibility of having a similar sensor placement and arrangement for all subjects while being able to accurately predict the biomechanical passivity index through the analysis of muscle networks. These results can be incorporated into a new family of controllers for pHRI systems, such as exoskeletons and rehabilitative robots, that rely heavily on the flow of energy from the robot to human biomechanics. The outcomes can address the current need explored in the literature for passivity-based stabilizers for these pHRI systems by reducing the conservatism of such closed-loop control platforms. This paper addresses the existing challenges of force-rendering transparency of such systems. Future work will include incorporating the proposed BioPI estimation technique into the design of existing passivity-based stabilizers for pHRI systems.