Morphological Computation of Haptic Perception of a Controllable Stiffness Probe

When people are asked to palpate a novel soft object to discern its physical properties such as texture, elasticity, and even non-homogeneity, they not only regulate probing behaviors, but also the co-contraction level of antagonistic muscles to control the mechanical impedance of fingers. It is suspected that such behavior tries to enhance haptic perception by regulating the function of mechanoreceptors at different depths of the fingertips and proprioceptive sensors such as tendon and spindle sensors located in muscles. In this paper, we designed and fabricated a novel two-degree of freedom variable stiffness indentation probe to investigate whether the regulation of internal stiffness, indentation, and probe sweeping velocity (PSV) variables affect the accuracy of the depth estimation of stiff inclusions in an artificial silicon phantom using information gain metrics. Our experimental results provide new insights into not only the biological phenomena of haptic perception but also new opportunities to design and control soft robotic probes.


Introduction
Physical examination of soft objects to identify hidden mechanical features can be seen in a variety of areas like minimally invasive surgery, medical physical examination, security, quality assurance in food industry, entertainment, etc. Manual examination involves variation of both behavioral and internal impedance levels of the fingers [1,2], because the interplay between motor control and internal mechanics (muscles and reflexes) [3,4] play an important role in both action and perception [5]. The notion-morphological computation-in soft robotics and biological systems views the mechanical circuits in the embodiment as a computational resource for both perception and action. Since the mechanical impedance of the embodiment changes the functionality of those mechanical circuits, it is important to understand its role in regulating perception and action of soft robots. The concept of internal impedance control in dynamic systems was first laid out by Hogan in [6]. However, its applications have been largely limited to control of action than understanding its role in perception or action-perception coupling in biological systems. For instance, impedance control principles have been widely applied in robotics for metastable walking [7], tele-operated excavation [8], safe interaction with humans [9,10], and to enhance stability of prosthetic devices [11]. They take advantage of passive body dynamics and the interaction with environment to achieve required goals without high-level cognitive processes [12]. Fig 1 illustrates how perception and action are coupled when they share a common embodiment. In biological systems, proprioceptive sensors such as spindle sensors (sense the amount and speed of muscle contraction) and tendon sensors (sense force) are located in the very muscles that are used to actuate joints. Therefore, controllers in the central nervous system perceive the environment depending on the actuation state of the muscles and muscle actuation in turn depends on perception [13]. Even in tactile perception [14], different types of mechanoreceptors are positioned at different depths in the dermis [15] to exploit different features of tissue dynamics. How this interplay among muscle actuation behavior, environment, and co-contraction of antagonistic muscles affect the accuracy of proprioception is not well quantified yet [16].
Recently, we have explored the role of internal impedance control for maximizing the information gain in robotic embodied perception [17]. The analytical results and experimental evidence have signified that a robotic manipulator with controllable internal impedance can maximize information gain (measured using transfer entropy) by searching for an optimal stiffness under a non-linear relationship between the entropy of sensory information and the impedance of physical embodiment of the robot. This paper addresses the question as to whether a robotic probe with variable impedance, indentation, and probe sweeping velocity (PSV) can achieve a better estimation accuracy of a given environmental condition (i.e. depth of a nodule embedded inside a silicon phantom).
In order to address this question, we designed and fabricated a new soft robotic probe with a controllable stiffness Mckibben type joint. The setup for probing experiments on a soft silicon phantom to estimate the depth of an embedded nodule is described in the following section. Then the dynamic simulations were carried out to explore the individual and collective role of internal impedance, indentation level, and PSV in the measured torque response. In the experimental section, we discuss experimental results derived from 5625 probing trials, where we explain our main contributions-1) Evidence to show that the internal stiffness of the soft probe plays a statistically significant role in the accuracy of nodule's depth estimation, 2) A Bayesian learning framework to combine internal stiffness, indentation level, and PSV that maximizes the information gain, and 3) Experimental results to show that proposed algorithm can achieve on average 99% and 96% accuracy in estimating the nodule's depth in both active and passive perception respectively.

Design and Analytical Model of Variable Stiffness Probe Design
The overall design of the probe used in this experiment is depicted in Fig 2(a). It is composed of two rigid links-tip link with length, l 1 = 80 mm, and base link with length, l 2 = 70 mmmade from ABS plastics. The joint connecting between these two links are coupled with a variable stiffness mechanism comprising of two identical linear ENTEX No.3552 stock springs (0.24 N/mm) from Advanex Europe Ltd, to represent how antagonistic muscles control the stiffness of biological fingers. Both springs are situated in spring chambers inside the base link and suspended between the pivot joint (at the connecting location, at which the relative angle to the vertical axis of the tip link is zero) and the movable anchor ring through a micro-filament thread (see Fig 2(b)). The stiffness of the joint can be regulated by changing the position of the anchor ring, r a , which is controlled by a 50 mm-linear actuator L12-50-210-06-I from Firgelli Technologies Inc. An ATI Nano17 F/T transducer is mounted at the top-end of the base link to measure the torque during the interaction with soft silicon phantom. This represents how the tendon sensor is located at the top end of a natural muscle. The probe's indentation level, i is controlled by a 30mm-linear actuator (L12-30-100-06-I from Firgelli Technologies Inc) mounted at the top of the base link. The total length of this mechanism when at rest is denoted by l o and has the initial value of 140mm. Hence the total rest length of this probe when the angle of pivot joint, q = 0, is 290mm. The probe structure is mounted on a flipped ANT130 XY-stage with built-in speed regulator from Aerotech Inc. as shown in Fig 2(c), which allows the planar movement in x-and y-direction at a constant speed. During the interaction with environment, e.g. palpation, the torque around the x-axis is Embodiment mediates both perception and action during the interaction between an agent and the environment. The system interacts with the environment through its embodiment. The internal impedance required for accurate perception through its embodied sensor can differ from that required to take appropriate action. Likewise an action taken with appropriate impedance could affect the quality of perception of the environment.  measured at the end of the base link of the probe using an ATI Nano17 Force/Torque (F/T) transducer.
Here, we used soft silicon phantoms with an embedded hard nodule as the samples for the haptic perception experiments. Silicon phantom is made from a soft clear silicon elastomer gel RTV27905 from Techsil. The given chemical substances (Part A and B) were mixed in 1:1 ratio according to specification to fabricate the silicon phantom. According to the studies in human biology, the typical tumor in soft tissue has roughly spherical shape [18] and at T1 stage can vary in size between 0.5 to 20mm [19]. Hence a spherical plastic bead with a diameter of 15mm was embedded inside each phantom at different depths (see Fig 2(d)). In this experiment, three silicon phantoms were used as samples where the nodules were embedded at the depth of 2, 4, and 8mm from the top surface of the phantom to the top of a nodule.

Variable Joint Stiffness
The exploded view of the variable stiffness joint is shown in Fig 2(b). At rest (the angular displacement, q = 0 and the position of anchor ring, r a = 0), the rest length of both springs are equal and denoted by r. Both springs can be extended from their rest by changing the position of the anchor ring, r a . The change of the angular displacement of the joint, q, due to the passive interaction with the environment, depending on the direction, causes both springs to be compressed and extended by Δr p , where Δr p = qR. R is the radius of the pivot joint at which the microfilament is attached to. Hence the change of the length of the springs can be computed as following: Since both springs are identical, we can compute the force contributed from each spring to the probe's joint as follows:f Hence, the torque generated from both springs due to the change of joint's angular displacement and the position of the anchor ring is: wheref s1;2? is the force produced from each spring perpendicular to the rotational axis.
Therefore, the total torque developed due to both springs can be computed from Eqs (3) to (5) as follows: and the stiffness at the joint, K s , is the derivatives of torque produced with respect to the angular displacement of the pivot joint, q, from Eq (6) From Eqs (6) and (7), we can simulate the torque, τ s , as well as the resulting joint stiffness, K s , as function of q and r a . The simulation results shown in The resulting joint torque due to the changes of both parameters, q and r a , is shown in Fig 3(a). Here we can see that the extension of both springs by increasing r a results in a change of the landscape of the relationship between τ s and q. Taking derivatives of the simulated joint torque with respect to the angular displacement of the pivot joint results in the stiffness profile of the joint shown in Fig 3(b). The stiffness of the joint becomes almost linear as the anchor ring approaches its origin. Since r a can be controlled through the linear actuator, for the rest of this paper, we represent the joint's stiffness level in term of the position of the anchor ring, r a mm.

System's Dynamics
The description of variables used in the system depicted in Fig 4 is shown in Table 1. The joint coupled between the tip and base link comprises of a variable stiffness element, which is represented by a variable spring-damper system. According to the derivation shown in Eq (7), the variable joint's stiffness, K s (r a , q), is therefore a function of the position of anchor ring, r a , and the angular displacement of the pivot joint, q.
Equation of Motion. The interaction dynamics of the system can be derived as follows: where τ j is the torque at pivot joint. The torque, τ = (τ f , τ j ) T can be resolved from the force components in y-and z-direction, f = (F y , F z ) > , at the probe's tip during the interaction with soft silicon phantom. τ f denotes the measured torque at the end of base link. The descriptions of the variables used here are explained in Table 1. Note that in this paper, the trigonometric functions are abbreviated as follow: s = sin(q) and c = cos(q). In order to simplify the dynamic equations of the system, the restoring force of the silicon phantom on the probe during the interaction can be modeled using a linear spring-damper system, where the stiffness of the silicon depending on different depth levels of a hard nodule embedded in the sample phantom can be represented by varying the system's spring constant.
Assume that: 1. At rest (no contact between probe and sample phantom) the probe has length of L = l o + l 1 + l 2 + i.
2. The base of the probe is fixed directly above the sample phantom at distance ρ.
3. The probe has stiffness K s as a function of r a and q.
4. The soft silicon phantom has uncertain stiffness k t with Gaussian distribution.
5. The restoring force from the soft silicon phantom is in both y and z-direction. 6. The friction between the tip and soft phantom's surface is negligible.
7. The deformation of soft silicon phantom has a uniform curvature [20].
When the probe comes in contact with the sample phantom, both the probe and the phantom deform according to their relative stiffness as shown in Fig 4. We denote the depth of phantom sample deformation in y-and z-directions by e y and e z respectively. The probe length (compressed prismatically), U, can be expressed as a function of q as Since the base of the probe is fixed, the constraint The probe consists of two rigid manipulator links with a variable stiffness joint. In the experiment, τ f is measured at the end of base link of the probe using ATI Nano17. The nodule embedded inside the phantom is shown in red-black spherical ball. The physical properties of phantom is described using spring-damper system. The variables used in the system is described in Table 1.
is maintained. Furthermore, we assume a uniform curvature deformation (the magnitude of deformation in both y-and z-directions are equal) of the soft silicon phantom. Therefore: By substituting Eq (9) in Eq (10), we obtain: The restoring force from the soft silicon phantom on the probe in both directions can be modeled as a spring-damper system as follows: Substituting Eqs (12)- (14) in Eqs (15) and (16), we obtain: Therefore the force component due to the interaction with soft silicon phantom at the tip is F y and F z are net force in y-and z-direction. f t is the translational force in probing direction. € y denotes the translational acceleration. Note that f y and f z contain variables, dependent on the terms indicated, since it is a function of the random variable k t . The terms r a ; q; _ q; and i, can be thought of as parameters to the distribution. Adjusting any of these may have an effect on the information in samples of f.
Torque measurement model. In the design of the probe shown in Fig 2(a), torque, τ f , is measured at the end of the base link. The Jacobian matrix, J of the system can therefore be expressed as: We get the equations of torque resulting from the interaction with the soft silicon phantom as follows: Hence the torque measured by the ATI Nano17 transducer mounted at the end of the base link can be derived as: Simulation According to Eqs (8) and (23), the torque response due to the interaction between the probe and soft tissue is dependent on the soft silicon phantom's stiffness k t , probe's stiffness r a , and indentation level i. Here we explore how different probing conditions such as: probe's joint stiffness, indentation, and PSV would affect the distribution of torque response at the probe's base during interaction with different phantom stiffness levels.
The expected value of the stiffness of the phantom, k m to , is identified to be 65 N/m with a standard deviation of 13.2 N/m. The source of uncertainty in the simulation is limited to that from the phantom's stiffness. Based on the previous studies [1,21], we can approximate the variability of the phantom's stiffness to be a Gaussian distribution, k t $ Nðk m t ; k s t 2 Þ, with expected value, k m t , and standard deviation, k s t . The changes in phantom's stiffness, Dk m t ¼ k m tn À k m to = 10, 20, 30, and 40 N/m, represent the presence of a hard nodule at different depths respectively. The length of the nodule can be viewed as the contact duration with the probe; hence the longer this is, the slower the PSV (v probe ). In the simulation, v probe is classified in three levels, namely: slow (10 mm/s), medium (20 mm/s), fast (30 mm/s).
The models of tissue's stiffness, in which the nodule is present, are given by where L t = 200 mm, l n = 15 mm, and t represent the length of the simulated phantom along the probing path, the diameter of the simulated nodule, and the simulation time respectively. t i and t f denote the time at which the probe's tip first contacts and leaves the nodule's area on the phantom's surface respectively. We simulated the dynamic torque response, τ f , during interaction between the probe and the phantom under different conditions specified in Table 2 As can be seen in Fig 5(b), τ f increases as the internal stiffness of the probe (controlled by the position of anchor ring, r a ) increases from r a = 0 to 4 mm. After that, τ f tends to settle down. The increase of the probe's indentation level also elevates the torque responses as shown in Fig  5(c). In Fig 5(d), the influence of the PSV, v probe , on the torque response, given the fixed values for the rest of the simulation parameters, cannot be visually assessed. Therefore, we have applied a statistical method to determine this. Since the simulated torque response is normally distributed (this was tested using Kolmogrov-Smirnov test for normality), we can implement ANOVA (Analysis of Variance) test. The result from the test signifies no statistically significant difference between these torque response distributions (p-value > 0.05). Therefore, the change in PSV, v probe does not significantly affect the torque response. From each torque response profile measured during palpation, we extracted the maximum torque at the location of simulated hard nodule as shown in 'circle' in These simulation results predict that the torque felt at the base of the probe can be controlled using probe stiffness, indentation level and PSV during the interaction with a soft tissue. The relationship between the torque felt and the combinations of probe's stiffness, indentation level, and PSV are non-linear. Furthermore, in reality the variability present in the system is non-deterministic and may arise from several sources such as the probing behavior, environment, and the probe itself. These raise the question as to how we can exploit these non-linear relationships to enhance the interpretation of the features in the environment using proprioceptive feedback from the torque sensor mounted at the base of the probe (representing how the tendon sensor is located in natural muscles). Since the relationship is stochastic and non- linear, the best way to preserve the interaction information is to present the relationship in the form of a probabilistic distribution. It provides us the opportunity to implement an appropriate stochastic machine learning technique to understand the role of individual factor in the interpretation of haptic perception.

Experiments and Results
In the experiment, we use the controllable stiffness robotic probe described earlier to derive deeper insights into the influence of the variation of combinations of probe's internal stiffness, indentation level, and PSV on the real-time estimation of the depth of a hard nodule embedded in soft silicon phantom. We explore whether a probe with controllable stiffness, indentation level, and PSV can exploit its past experience of palpation by varying its own internal stiffness, indentation level, and PSV to estimate the depth of embedded nodule inside a soft silicon phantom.

Palpation Memory Primitives
The palpation experiences during interaction with the soft silicon phantoms with a nodule embedded at different depths can be presented in the form of a probabilistic representation, which hereafter is referred to as 'memory primitives'. The memory primitives were built from torque measurements, τ f , for different r a , v probe and i, over multiple palpation learning trials.
The probe mounted under the XY-table was programmed to palpate in a straight line along the probing path over the soft silicon phantom's exposed surface shown in Fig 2(d). During each palpation trial, the torque response due to the interaction is measured at the base around the x-axis at 1000 Hz. Readings of r a , i, v probe , and τ f were recorded using LabVIEW2012 software, National Instruments Corp, through the data acquisition cards. Data processing was carried out using MatLAB R2013b, The MathWorks, Inc.
In order to construct the primitives for this experiment, we conducted palpation experiments across 5 probe stiffness levels, r a , 5 indentation levels, i, 3 levels of PSV, v probe , and 3 levels of nodule depths, d. The experimental conditions are shown in Table 3. These amount to 225 unique interaction conditions. For each given combination, 25 palpation trials were repeated to allow the formation of distribution of peak torque response.
Each measured torque signal from the F/T transducer is de-noised for 5 levels using wavelet decomposition technique with a Daubechie's db10 mother wavelet. The peak torque at the nodule's location is then extracted from each de-noised signal. The peak torque distribution given different combinations of the probe's stiffness, indentation level, and the probing speed for different nodule's depth level, P(τ f |d, r a , i, v probe ), can be constructed from all 25 trials by fitting a normal distribution to the data. Here, only 81 interaction conditions are depicted as examples of memory primitives as shown in Fig 6(a), 6(b) and 6(c). Experimental results and analysis The non-linear relationship between the measured torque probability distribution and different combinations of probe's stiffness, indentation level, and PSV can be used in an appropriate machine learning algorithms to enhance the accuracy in nodule's depth estimation. In this paper, interpretation of real-time torque measurements during palpation was done using memory primitives (conditional probability density functions) in a Bayesian Inferencing framework as given by P t ðdjt f Þ ¼ Pðt f jd; r a ; i; v probe ÞP tÀ1 ðdÞ P m n¼1 Pðt f jd n ; r a ; i; v probe ÞP tÀ1 ðd n Þ ; ð25Þ where P t (d|τ f ) is the posterior probability of nodule's depth given τ f . t denotes the current estimation iteration. P t−1 (d) represents the prior distribution of d. P(τ f |d, r a , i, v probe ) refers to the likelihood probability distribution of torque, given d, r a , i, and v probe . n indicates the index of d. And m = 3 is the number of possible classes for nodule depth. Bayesian Inference in Estimating the Nodule's depth. Table 3 shows different combinations of probe-soft tissue interaction conditions used to assess the performance of the Bayesian Inferencing algorithm 1 to estimate d across 5 iterations. The assessment of the nodule's depth estimation was repeated for 10 trials in order to obtain the average accuracy of estimation. In each assessment trial, the memory primitives given different combinations of probe's interacting conditions were constructed from 25 randomly chosen learning trials. for each iteration t 2 1..5 do 5 Retrieve and process new τ f t from the sensor reading, given known probing bahavior { r a , i, v probe }.; Compute P(τ f t |d, r a , i, v probe ) from φ.; 7 Recall prior distribution of hypothesis of nodule's depth P t−1 (d).; 8 Compute P t (d|τ f ) using Eq 25.; 9 Store posterior distribution as a prior distribution for the next iteration.; The overall accuracy across 10 assessment trials in nodule depth estimation after each iteration and those for different nodule depth levels are shown in Fig 8. On average, the overall nodule depth estimation accuracy increases from approximately 91% with standard deviation of 3.23% at the first iteration (t = 1) to 96% with standard deviation of 1.8% at (t = 5). We witness that the estimation accuracy of nodule depth decreases as the nodule is buried deeper from the exposed surface from 99.3% at d = 2 mm to 95.2% and 94.2% at d = 4 and 8 mm, respectively. Higher iteration numbers cause the expected values of estimation accuracy for all depth ranges to increase and the standard deviations to decrease.
Further statistical analysis was performed to investigate the significance of r a , i, and v probe , on the depth estimation accuracy across all assessment trials. The Kolmogrov-Smirnov test showed that the nodule's depth estimation was not normally distributed p > 0.05. Hence, the conventional Analysis of Variance (ANOVA) could not be performed. Therefore, a nonparametric Kruskal-Wallis method was applied. The resulting p-values from the test were 0.0002, 0.4715, and 0.7394 for probe's stiffness, r a , indentation, i, and PSV, v probe , respectively over the course of 10 assessment trials across different combinations. Therefore, probe's stiffness, r a , has statistically significant contribution towards the nodule's depth estimation accuracy (p-value <0.05); while the variation of i and v probe do not have a significant influence.
The average accuracy in the estimation of nodule's depth given different probe's stiffness, r a , indentation level, i, and PSV, v probe across 10 assessment trials are shown in Tables 4, 5 and 6, respectively. We can see that the average nodule's depth estimation is slightly less accurate for  the case when the joint of the probe is both stiffest and most relaxed; while it is most accurate when r a = 12 mm. This suggests that the regulation of body's stiffness matters when making an estimation about the environment. Nonetheless, it is important to notice that the efficacy of r a = 12 mm varies depending on the probing speed and the indentation used. Furthermore, Fig 7 suggests that the posterior distribution of nodule's depth estimation can converge at different rates depending not only on the combination of r a , i, and v probe , but also on the noise level of real-time sensor measurements in each iteration and the chosen memory primitives in the likelihood function. This leads to the question as to whether we can determine the sufficiency of the number of iteration/exploration required to make an accurate estimation of the nodule depth. We can address this question by computing the information transfer entropy in each iteration. This will be addressed in detail in the next section.
Kullback-Liebler Transfer Entropy. In general, the common influences of multiple coupled systems and factors can be quantified through the directed information exchanges by measuring the information transfer entropy, also known as relative entropy [22]. For example, we can assign the combination of internal stiffness, indentation level, and PSV; and the torque sensor reading to be random variables (RV-A) and (RV-B) respectively. While the mutual information of two coupled variables between RV-A and RV-B does not change with the exchanges of variables; the transfer entropy from RV-A to RV-B is not identical to that from RV-B to RV-A. Transfer entropy can be quantified using Kullback-Liebler (KL) divergence.
To be more specific, KL-divergence can be used to assess whether further information regarding the nodule's depth estimation can be gained by taking another action (further iteration in Bayesian nodule's depth estimation procedure). If we consider a set of P t (d|τ f ) computed at the end of each Bayesian iteration as the hypothesis of the depth estimation, its entropy for a given torque measurement, τ f , is dependent on a set of probe's stiffness, indentation level, and PSV, {r a , i, and v probe }. KL-divergence defined in Eq (26) represents the additional information gained, G, about the relationship between the hypothesis of depth estimation, P t (d), and τ f across iterations of Bayesian Inference as well as across different sets the probe's stiffness, indentation level, and PSV. Therefore, KL-divergence is a good measure to quantify the gain of different actions underlying the changes in the behavior.
P t (d|τ f ) represents the probability distribution of depth estimation which is obtained from the Bayesian inference shown in Eq (25) at t th iteration, and P t = 0 (d) represents the base hypothesis about the nodule's depth estimation. Bayesian Inference in Estimating the Nodule's depth with Kullback-Liebler Divergence. In addition to the Bayesian Inference method for estimating the depth of the nodule from the real-time captured torque data, here KL-divergence is implemented at the end of each Bayesian iteration to determine whether further measurement is required to accurately estimate the nodule's depth. This additional process is carried out by computing the correlation distance, δ, between information gain from the current hypothesis, G t , and that from the prior hypothesis, G t−1 , in relation to the base prior distribution, P t = 0 (d). No further measurement is necessary to make an accurate estimation when the correlation distance, δ, is less than the empirically specified threshold of T = 0.0005. This signifies that there is negligible change in the information gained across the iterations and the distribution of the nodule's depth estimation hypothesis has converged. The depth estimation procedure is shown in Algorithm 2. Similar to the process presented in previous section, the assessment of the nodule's depth estimation is repeated for 10 assessment trials to obtain the average accuracy in the estimation. In each assessment trial, the memory primitives given different combinations of probe's interacting conditions are constructed from 25 randomly chosen learning trials.
Algorithm 2: Nodule's depth estimation algorithm using Bayesian Inference and KL divergence 1 function DepthEstimationKL(τ f (d r , r a , i, v probe )); Input: Real time torque reading, τ f (d r , r a , i, v probe ) Output: Depth estimation accuracy 2 Define P t = 0 (d) as a flat distribution across different d; 3 for each combination of { r a , i, v probe }, and actual nodule's depth, d r do 4 δ = 1 Initialize correlation distance to 1; 5 t = 0 Initialize the number of iteration to 0; 6 G 0 = 0; Initialize the information gain at t = 0 to 0; 7 while δ > T do 8 t = t + 1; 9 Follow Step 5-9 in Algorithm 1; 10 Compute G t using Eq (26); 11 Compute δ between G t and G t−1 . With the implementation of KL-divergence in addition to the Bayesian Inference algorithm, the nodule's depth estimation process requires on average of only 2.8 iterations with standard deviation of 1.2 iterations to converge. While the number of iterations required for convergence is kept to minimum; the nodule's depth estimation accuracy still reaches within the comparable range to that with fixed 5-iterations in the inferencing algorithm presented in Algorithm 1. On average the overall depth estimation accuracy is approximately 96.2% as shown in Fig 9 (orange bars). The accuracy of nodule's depth estimation for each actual depth are approximately 98.4%, 95.3%, and 94% for d r = 2, 4, and 8 mm respectively. These results show that this method minimizes the number of explorations needed to make an accurate estimation about the depth of the nodule. Therefore, we can conclude that Bayesian Inference together with KL-divergence provides a real-time framework to estimate the convergence to an optimal estimate of nodule depth in the sense of information gain.
Bayesian Inference in Estimating the Nodule's depth with Dynamic Probing. So far, the experiments involved keeping r a , i, and v probe constant given a set of probing iterations in a trial. However, biological counterparts like humans regulate the internal impedance more like a random variable within a given probing attempt. Therefore, we conducted further experiments to explore whether we can enhance the nodule's depth estimation accuracy by allowing changes in the combination of probe's stiffness, indentation level, and PSV across trials. This allows the estimation process to explore in multiple search spaces (memory primitives).
The nodule's depth estimation result from the implementation of the Bayesian Inference with dynamic probing shown in Algorithm 3 are shown in blue bar in Fig 9. The overall average accuracy from 100 trials using this algorithm as the estimation hypothesis converges reaches 99% with standard deviation of 0.5%. For each individual depth level, the average estimation accuracy were all significantly higher; while the corresponding standard deviations were lower compared to those from the Bayesian Inference with static set of probe's stiffness, indentation level, and PSV combinations. The result from this assessment also confirms our initial hypothesis that the average number of iterations required to perform accurate depth estimation is kept to a minimum at approximately 3 iterations with standard deviation of 1.3 iterations.
A preliminary experiment was carried out with 1 human subject in the same probing task. In order to have comparative basis between human and robotic experiments, the subject was blindfolded and asked to palpate the same set of soft phantoms to estimate the depth of the embedded hard nodule. The muscle activity caused by the stiffness regulation of the finger was captured using electromyography (EMG) signal at Flexor digitorum superficialis (FDS) and Extensor digitorum communis (EDC) [23]. The example of both normalised EMG signals are shown in Fig 10(a) and 10(b). The combination of both EMG signals results in the co-contraction activity of the muscle pair. The muscle co-contraction pattern of a human subject during manual palpation is shown in Fig 10(c). The peaks (circled in red in Fig 10(c)) indicate the coactivation of the muscle pair, when both FDS and EDC are contracted. Similar to the robotic experiment with the implementation of Algorithm 3, the human's co-contraction pattern also demonstrates variability. The variability of both EMG signals arises from the changes in the finger's stiffness level of which the activation of FDS and EDC muscle pair are responsible for. The regulation of internal stiffness level of the human's finger can be correlated with the change in the joint's stiffness of the robotic probe. Our results based on robotic experiments explain the reason behind this activity.

Discussion and Conclusion
This paper investigated both individual and collective role of the probe's internal stiffness, indentation level, and PSV in the accuracy in interpreting and estimating an environmental feature (depth of a nodule) by controlling a soft probe. The soft robotic probe comprised of a variable stiffness joint and an indentation level control mechanism. The probe structure was mounted under an XY-stage allowing the planar movement. Firstly, we simulated the dynamics of palpation using the designed probe on a simulated soft silicon phantom to observe the interaction between them under different combinations of probe's internal stiffness, indentation level, and PSV.
The results from the simulation suggest that 1) the torque felt at the base of the probe can be controlled using different combinations of probe's stiffness, indentation level, and probing speed and 2) the relationships between the torque measured, the stiffness of the soft silicon phantom, and the combination of probe's internal stiffness, indentation level, and PSV are non-linear. While a variability in the simulated system from the Gaussian distribution of the phantom's stiffness is pre-defined; the variability in such a system in reality is non-deterministic and can arise from multiple sources. These brought into a question as to how we can use these non-linear relationships in the experiment to enhance the estimation of the environmental features.
In the experiment, we investigated the question as to how the probe with controllable stiffness, indentation level, and PSV can exploit its past experience of palpation to estimate the depth of a nodule embedded inside a soft tissue in real time. The non-linear relationship between the probe's measured torque, its internal variables, and the environment (depth of nodule in silicon phantom) were presented in the form of a probabilistic distribution given different combinations of probe's internal stiffness, indentation level, and PSV. In this paper, we referred to these conditional probability distributions as 'memory primitives'. These 'memory primitives' functioned as likelihood functions in a Bayesian framework to estimate the depth of a nodule in the soft tissue phantom. The memory primitives were constructed from three levels of PSV, five levels of indentation, and five levels of joint stiffness, for three nodule's depth levels. In total 5625 probing trials were performed using this automated experimental setup.
In conclusion, the implementation of Bayesian Inference allows the algorithm to accurately estimate the depth of a nodule from the measured torque real-time. Furthermore, KL-divergence was introduced to determine whether further iteration of measurement is required to make an accurate estimation by comparing the information gained in the current iteration to that of the previous iteration. It was shown that on average the estimation processes using Algorithm 2 and 3 require approximately 3 iterations to converge in order to obtain comparable and better (in the latter) estimation accuracy. Finally, allowing the combination of probe's internal stiffness, indentation level, and PSV to randomly vary across iterations (allowing exploration in multiple memory primitives in each nodule's depth estimation process), resulted in a convergence to the global optimum with a minimum number of iterations. We showed that, this could enhance the average depth estimation accuracy to almost 100% with higher repeatability (smaller standard deviation).
The insights from this study sheds light on the practice of manual and robotic assisted palpation of soft tissue to locate T-1 stage tumors in biological tissues [18]. Medical literature shows that T-1 stage tumors can be modeled as spherical shape hard nodules. Since the focus of this paper is to understand the importance of the internal impedance of the probe in detecting a hard nodule in a soft tissue, we limited the study to a spherical acrylic hard nodule buried at depths upto 8mm. This scenario represents the conditions of a typical manual tumor localizing procedure for a T-1 stage tumor. Even within this range of depths, Fig 9 shows how the accuracy of nodule's depth estimation decreases as the depth of the nodule increases. Nodules located deeper in the tissue would require higher indentation forces potentially causing damage to the tissue. However, future studies could be done with different shapes and materials of nodules buried deeper in the tissue. This paper has provided important explanations about the role of morphological computation in haptic based probing of a soft object, as well as providing guidelines to design and control variable stiffness probes for physical examination. Certainly, the operational implementation of this probe should be further developed depending on different applications. However, the fact that controllable internal stiffness helps to gain proprioception information is still valid in such a tool. However, additional complexities arising from factors such as variable friction and irregular surface conditions not addressed in this paper should be further examined. Future studies will also involve temporal control of probe stiffness, indentation, and speed to better understand diverse probing strategies used by different classes of human participants as seen in [1].