Fused ultrasound and electromyography-driven neuromuscular model to improve plantarflexion moment prediction across walking speeds

Background Improving the prediction ability of a human-machine interface (HMI) is critical to accomplish a bio-inspired or model-based control strategy for rehabilitation interventions, which are of increased interest to assist limb function post neurological injuries. A fundamental role of the HMI is to accurately predict human intent by mapping signals from a mechanical sensor or surface electromyography (sEMG) sensor. These sensors are limited to measuring the resulting limb force or movement or the neural signal evoking the force. As the intermediate mapping in the HMI also depends on muscle contractility, a motivation exists to include architectural features of the muscle as surrogates of dynamic muscle movement, thus further improving the HMI’s prediction accuracy. Objective The purpose of this study is to investigate a non-invasive sEMG and ultrasound (US) imaging-driven Hill-type neuromuscular model (HNM) for net ankle joint plantarflexion moment prediction. We hypothesize that the fusion of signals from sEMG and US imaging results in a more accurate net plantarflexion moment prediction than sole sEMG or US imaging. Methods Ten young non-disabled participants walked on a treadmill at speeds of 0.50, 0.75, 1.00, 1.25, and 1.50 m/s. The proposed HNM consists of two muscle-tendon units. The muscle activation for each unit was calculated as a weighted summation of the normalized sEMG signal and normalized muscle thickness signal from US imaging. The HNM calibration was performed under both single-speed mode and inter-speed mode, and then the calibrated HNM was validated across all walking speeds. Results On average, the normalized moment prediction root mean square error was reduced by 14.58 % (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=0.012$$\end{document}p=0.012) and 36.79 % (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<0.001$$\end{document}p<0.001) with the proposed HNM when compared to sEMG-driven and US imaging-driven HNMs, respectively. Also, the calibrated models with data from the inter-speed mode were more robust than those from single-speed modes for the moment prediction. Conclusions The proposed sEMG-US imaging-driven HNM can significantly improve the net plantarflexion moment prediction accuracy across multiple walking speeds. The findings imply that the proposed HNM can be potentially used in bio-inspired control strategies for rehabilitative devices due to its superior prediction. Supplementary Information The online version contains supplementary material available at 10.1186/s12984-022-01061-z.

Page 2 of 20 Zhang et al. Journal of NeuroEngineering and Rehabilitation (2022) 19:86 Background Locomotion mobility accounts for a dominant part of human activities of daily living, like moving around the home and community, going to work or school, doing errands, visiting friends, etc. The human lower extremity plays an essential role in achieving locomotion mobility. The human ankle plantarflexors generate a large burst of "push-off " mechanical power during the late stance phase of walking, enabling forward and upward acceleration of the body's center of mass. Due to neurological disorders or injuries like spinal cord injury, stroke, and multiple sclerosis, the weakened function or dysfunction of plantarflexors is likely to cause a dramatic decrease in the "push-off " power. Consequently, these mobility disorders impair walking function and cause poor energy economy [1], as well as disrupt both physical and emotional well-being [2]. Modern neurorehabilitation devices, such as powered ankle exoskeletons [3][4][5], soft exosuits [6,7], and functional electrical stimulation [8][9][10][11][12], may use assist-asneeded control to actively engage and maximize recovery of users with mobility impairments [13,14]. In turn, the efficacy of the control strategy depends on the accurate determination of continuous human volitional movement intent (net joint moment). Mechanical sensors, like force or torque sensors, installed on a rigid and bulky frame, have been often used to measure the human intent for joints without direct interaction to the ground, but limit the system's wearability. Also, inaccuracies may creep in easily due to the inevitable misalignment between the bionic joint center and human joint center, which may introduce undesired interaction force [15,16]. Usually, it is challenging to directly measure the net ankle joint moment during walking overground with conventional force or torque sensors setup. The standard way to measure the net ankle joint moment uses a motion capture system, ground reaction force (GRF), and inverse dynamics (ID) calculations. However, there are two shortcomings of the standard approach. First, the setup is constrained to a lab environment, and not applicable for field testing. Second, the results from ID do not reveal how skeletal muscles perform during walking from a neuromuscular perspective. Therefore, a forward dynamics approach based on a neuromuscular model would be very instrumental when a motion capture system and GRF data are unavailable.
Surface electromyography (sEMG) measures electrical potentials during asynchronous muscle neurons firings, and its amplitude and frequency positively relate to muscle activation levels. Therefore, sEMG-derived signals can be used in a Hill-type neuromuscular model (HNM) or to train a machine learning approach (model-free) to predict volitional joint moment [17,18] and angular position [19][20][21]. However, sEMG signals suffer from interference or cross-talking from the adjacent muscles, and an inability of measuring activations of deep-layer muscles [22,23]. Alternatively, two-dimensional brightness mode (B-mode) ultrasound (US) imaging allows one to see the musculature of the targeted muscle in vivo. Due to its ability to directly visualize superficial and deeplayer muscles, US imaging may work as an alternative methodology to predict joint motion or motion intent. Potentially, US imaging overcomes the shortcomings of the sEMG measurements. Most frequently used architectural features from US images include pennation angle (PA) [24,25], fascicle length (FL) [26,27], muscle thickness (MT) [28,29], and cross-sectional area [30]. These features have been correlated with the joint kinetics and kinematics during isometric or isokinetic joint motion by using HNM-based or model-free approaches [25,[31][32][33].
Motivation also exists to use a dual-modal approach that combines the measurement from sEMG and US imaging. Potentially sEMG signals and US imaging provide complementary information, and the combination between them may (1) mitigate any cross-talking effect from neighboring sEMG signals and (2) lower US imaging-derived features' drift due to the accumulated errors from cyclic joint movement. Our recent studies have shown the advantages of using the dual-modal approach over uni-modal bio-signals (sEMG or US imaging) for ankle joint moment/motion prediction under isometric/ dynamic dorsiflexion studies [21,25,27] and isometric plantarflexion study [34,35]. Similarly, Dick et al. used both sEMG and US imaging to predict plantarflexion force during dynamic cycling tasks [32]. However, the use of the dual-modal bio-signals during more complex functional tasks, like walking across different speeds remains unexplored.
In this study, we investigated a dual-modal approach that takes both processed sEMG and MT from US imaging as inputs to a modified HNM, named sEMG-US imaging-driven HNM, to predict net ankle joint plantarflexion moment during the walking stance phase across multiple speeds. We hypothesize that (1) the proposed HNM will achieve a better net plantarflexion moment prediction accuracy than sEMG-driven and Page 3 of 20 Zhang et al. Journal of NeuroEngineering and Rehabilitation (2022) 19:86 US imaging-driven HNMs, (2) the net plantarflexion moment prediction performance is more robust if data collected from multiple speeds are included in HNMs' calibration procedure.
In previous US imaging studies, visualized architectural features, such as PA and FL, require high-resolution US imaging, which can be significantly affected by a US transducer placement site on the muscle. To mitigate the requirement of high-resolution US imaging, US imagingderived signals such as echogenicity [27,36,37], tissue displacement [38], and MT [29], are more preferable to correlate with muscle or joint mechanical functions. Because ankle plantarflexors: lateral and medial gastrocnemius and soleus muscles (LGS, MGS, and SOL) are not accessible in the same US image plane, this study chose to focus on LGS and SOL, which are in the same plane, and tracked their MT change during the walking experiments across multiple speeds. Therefore, the contributions of the paper are: (1) the use of US imaging-derived MT changes and sEMG-derived changes of LGS and SOL muscles as surrogate variables of muscle activations during walking, (2) development of a modified HNM that uses both sEMG and US imaging as inputs for singlespeed modes and an inter-speed mode, and (3) evaluation of the sEMG-US imaging-driven HNM's robustness across multiple walking speeds.

Hill-type neuromuscular model of ankle joint
Below, we propose a modified sEMG-US imaging-driven HNM to directly build a relationship between the joint net PF moment and sEMG-US imaging-derived surrogate signals of both LGS and SOL muscles. There are three sub-models involved in the HNM: (a) sEMG-US imaging-derived weighted muscle activation model, (b) muscle-tendon unit geometry model, and (c) muscle contraction dynamic model.

Weighted muscle activation model
The neural activation at t k time instant N i (t k ), i = 1, 2 , k = 1, 2, 3, ... , for SOL and LGS muscles, respectively, considers the electromechanical delay (EMD), τ , between the onset of an sEMG signal and a muscle contraction and utilizes a second-order recursive filter and is defined as [39] where EMD, τ , is usually between 30 ms and 120 ms. u i (t k ) represents the sEMG's linear envelope normalized to the peak value of the specific task (with a constant peak value across all walking speeds on each subject in this study). The sEMG's linear envelope was derived after raw sEMG signals' band-pass filtering, full-wave rectification, and low-pass filtering. α i , β 1i , and β 2i are coefficients that define the recursive filter's dynamics of each muscle [39], and the following set of constrains is employed to reach a positive stable solution, i.e.
A nonlinear relationship between neural activation N i (t k ) and corresponding muscle activation a 1i (t k ) is given as [39] where A i represents the nonlinear shape factor for each muscle, which is allowed to vary between − 3 and 0, with A i = −3 being a nonlinear and A i = 0 being a linear relationship.
According to the reported results in [40][41][42][43], the MT change was found to correlate with the muscle contraction level or muscle activation through a linear function or piece-wise linear function. Therefore, in this work, the second part of the LGS or SOL muscle activation is calculated from the US imaging-derived MT change. The MT values of LGS and SOL muscles are denoted as MT i (t k ) . According to the preliminary results in [44], there is a positive relationship between MT change and the targeted muscle contraction level during the walking stance phase. After taking normalization of the MT with respect to the lower bound (at rest state) and upper bound (peak value of the specific task) on each participant, the US imaging-derived muscle activation, a 2i (t k ) , is proposed as where MT i max and MT i min represent the constant subjective MT values when the LGS and SOL muscles are at the walking task-specific maximum voluntary contraction and complete rest condition, respectively. For each individual, MT i max and MT i min are set as consistent values across different walking speeds. This normalization guarantees that the US imaging-derived muscle activations vary between 0 and 1.
By introducing an allocation gain between the sEMGand US imaging-derived muscle activations, the synthesized/weighted muscle activation levels for LGS or SOL muscle, a i (t k ) , can be represented as where δ i ∈ [0, 1] represents the muscle activation allocation gain for LGS ( i = 1 ) and SOL ( i = 2 ) muscles. (2)

Muscle-tendon unit geometry model
As presented in Fig. 1b, consider the ankle joint rotation center in the sagittal plane as O, the proximal and distal osteotendinous junction points of each each muscletendon unit (MTU) near the knee joint and heel as A i and B i , and the angle between OA i and OB i as q(t k ) , then each MTU length, l mt i (t k ) , is represented as the distance between A i and B i , and is calculated as where l OA i and l OB i represent the distance of OA i and OB i obtained from OpenSim (National Institutes of Health for Biomedical Computation, Stanford, USA) [45], respectively. A generic OpenSim model (gait2392) was linearly scaled to each participant in OpenSim version 4.1 [45] per [46]. According to the law of sines, the moment arm of each MTU, r mt i (t k ) , is calculated as The MTU generates contraction force only when it is stretched, which indicates the current tendon length l t i (t k ) is equal to or longer than the tendon slack length l sk t i . From the perspective of muscle contraction dynamics that is shown in Fig. 1c, the overall MTU length, l mt i (t k ) , could also be expressed as where l t i (t k ) and l m i (t k ) represent the current tendon length and muscle fascicle length for LGS and SOL muscles, respectively. φ i (t k ) represents the pennation angle that changes with instantaneous l m i (t k ) . By assuming the individual muscle belly has a relatively small MT change and volume change [47,48], φ i (t k ) can be approximately calculated as where the term φ 0 i denotes the constant pennation angle when the muscle is at optimal fascicle length l 0 m i . The work in [49,50] has shown that the optimal fascicle length increases as muscle activation decreases. Due to this coupling, the following relationship between the muscle activation and corresponding optimal fascicle length is given as where is the rate of change in the optimal fascicle length, and it is selected as 0.15 [39]. l 0 mo i is the optimal fascicle length at the walking task-specific maximum voluntary contraction, and l 0 m i (t k ) is the optimal fascicle length at time t k and muscle activation a i (t k ).
From the preliminary results of plantarflexors' US imaging in [44], it is very challenging to capture the entire fascicle length l m i (t k ) of LGS or SOL muscles by using the current US transducer with a small dimension (width of 38 mm). Therefore, the fascicle length parameters of either LGS or SOL muscles were not measured directly from the US images. Instead, an indirect calculation method was applied as detailed below. According to [51], the nonlinear nominal tendon force-tendon strain relationship is given as represents the tendon strain and F max i represents the muscle contraction force at the walking task-specific maximum voluntary contraction. By can be numerically computed based on the Runga-Kutta-Fehlberg algorithm. By substituting (6) and (9) to (8), the muscle fascicle length l m i (t k ) will be calculated, as well as the muscle fascicle velocity v m i (t k ) by taking the time derivative of l m i (t k ).

Muscle dynamic contraction model
In the sEMG-US imaging-driven HNM, the individual moment component produced by each muscle, M i (t k ) , can be represented as where r mt i (t k ) is defined in the above geometry model and F mt i (t k ) denotes the corresponding contraction force generated on each MTU and is represented as where φ i (t k ) represents the pennation angle defined above. F ce i (t k ) and F pe i (t k ) denote the corresponded forces generated by the parallelly located contractile element and the passive element, and can be calculated as In (14), F max i of the LGS and SOL muscles will be identified based on optimization algorithm in the HNM calibration procedures.
) , and f p i (l m i (t k )) denote the generic muscle contractile forcefascicle length, force-fascicle velocity, and passive elastic force-fascicle velocity curves. These curves were normalized to F max i , optimal fascicle length l 0 m i , and maximum fascicle contraction velocity v max m i . The explicit expres- ) , and f p i (l m i (t k )) can be found in [25,32,52].

Experimental protocol, data collection and pre-processing
The study was approved by the Institutional Review Board (IRB) at North Carolina State University (IRB number: 20602). Ten young participants (7M/3F, age: 25.4 ± 3.1 years, height: 1.77 ± 0.10 m, mass: 78.0 ± 21.1 kg) without any neuromuscular or orthopedic disorders, were recruited in this study. Every participant got familiar with the experimental procedures and signed an informed consent form before participation. Figure 2a summarizes the experimental setup for the participants to perform static anatomical poses and dynamic gait trials (walking speeds changing from 0.50 to 1.50 m/s), and Fig. 2b presents the workflow of the data processing and sEMG-US imaging-driven HNM calibration procedures. During all walking trials, threedimensional coordinates of 39 retro-reflective markers positioned on the participant's lower extremities following the instructions in [53] were recorded using a 12-camera motion capture system (Vicon Motion Systems Ltd, Los Angeles, CA, USA) at 100 Hz. The GRF signals were collected at 1000 Hz synchronously with makers trajectories using in-ground force plates (AMTI, Watertown, MA, USA) mounted on an instrumented treadmill (Bertec Corp., Columbus, OH, USA) through the commercial real-time data capture software Nexus 2.9. Both GRF signals and markers trajectories were low-pass filtered with a fourth-order Butterworth filter in Visual 3D software (C-Motion, Rockville, MD, USA), and the cut-off frequencies were set as 25 Hz and 6 Hz, respectively. The markers trajectories and GRF signals from static poses and dynamic gait trials were used for joints kinematics and net joint moment calculation by using inverse kinematics and ID algorithms in Visual 3D. During all trials, four wired sEMG sensors (SX230, Biometrics Ltd, Newport, UK) were attached to the shank skin through a double-sided tape to non-invasively record the electrical signals from tibialis anterior (TA), LGS, MGS, and SOL muscles at 1000 Hz synchronously through Nexus 2.9. The locations for sEMG sensors were determined following the instructions in [54]. A US transducer (38-mm of length, 6.4 MHz center frequency, L7.5SC Prodigy Probe, S-Sharp, Taiwan) setup was applied in the walking experiments as described in [34]. The US radio frequency (RF) data were recorded at 1000 frames per second synchronously with markers trajectories, GRF, and sEMG signals, using a pulse sequence trigger signal generated from Nexus 2.9. sEMG signals were band-pass filtered with the bandwidth between 20 Hz and 450 Hz, and then full-wave rectified and low-pass filtered with a cut-off frequency of 6 Hz. The resulting linear envelopes were normalized to the peak processed sEMG values obtained from all sets of dynamic gait trials with different speeds. The US radio frequency data were converted to B-mode images through beamforming and logarithmic compression, and then a commercial US imaging processing toolbox UltraTrack [55] was applied to determine the temporal sequences of the LGS and SOL muscles' thickness based on the adaptive optical flow tracking algorithms. First, one region of interest (width by height: 336 × 400 pixels) that encompassed the LGS and SOL muscles was selected as the area between the superficial and deep aponeuroses. Then, 10 vertical lines were manually defined with evenly distributed distances in the LGS's region and SOL's region on the first US imaging frame for each recorded waking trial, as shown in Fig. 2b. Key-frame correction [55] was applied to minimize the time-related drift of MT's cyclical pattern over multiple gait cycles, where the key frames were selected to be at heel-strike and toe-off time points. After each correction, the new key frames' vertical lines positions were determined by applying an affine transformation to the key-frame before it. Finally, the temporal sequence of the mean value of those 10 vertical lines' lengths from LGS and SOL muscles were calculated for LGS's MT and SOL's MT, respectively. The MT signals were then low-pass filtered with a cut-off frequency of 30 Hz.
Before data collection, participants practiced dynamic walking steadily with all sensing devices attached to the lower extremities on the treadmill for at least 30 s for each speed. The duration of each walking trial was set as 2 min, and data from the middle 20 s were collected for analysis. In the current study, there were five different walking speeds, 0.50 m/s, 0.75 m/s, 1.00 m/s, 1.25 m/s, and 1.50 m/s, and the order was selected randomly for each participant. The participants were provided with at least 2 min for rest between two successive dynamic walking trials to avoid muscle fatigue.

Data post-processing HNM model calibration
To determine the values for a set of parameters, including the shape factor A i , the tendon slack length l sk t i , the muscle contraction force at the walking task-specific , and the muscle activation allocation gain δ i , the HNM model calibration was performed with the initial parameters obtained from the literature [56]. The summarized calibration process of the sEMG-US imaging-driven HNM is presented in Fig. 2b. Indeed, the model calibration is a nonlinear optimization problem, where the above parameters set needs to be found for minimizing the designated objective function. During calibration, the parameters were subsequently adjusted within the predefined boundaries, which ensured the muscle-tendon unit always operated within the physiological range [39]. Except for the allocation gains δ i were set the range between 0 and 1, other parameters were set with the lower and upper bounds as 50% and 150% of the literature-referred values, while the initial guess of the set of parameters were set as the middle value of each parameter. Other HNM parameters, like muscle-tendon unit length l mt i , optimal muscle fascicle length l 0 mo i , and optimal pennation angle φ 0 i were assigned using OpenSim and the scaling method as references. The Matlab function 'lsqcurvefit' with the Levenberg-Marquardt algorithm was applied to solve the nonlinear least-squares optimization problem until the following objective function was minimized during the calibration represents the net PF moment estimation from both LGS and SOL muscles during stance phase at time instant j, M w (j) represents the net PF moment measurement from ID at time instant j, and N represents the length of the data used for the calibration.
It should be noted that the HNM calibration is subjective, and for each participant, the calibration procedure was repeated for both single-speed modes and interspeed mode. In the current study, single-speed modes are defined as the HNM calibration only with data collected from one single speed (five single-speed modes here), while the inter-speed mode is defined as the HNM calibration with all data collected from five different speeds (one inter-speed mode here). In HNM calibration, collected data, including ankle joint's kinematics, kinetics, sEMG signals, and US imaging from the stance phase of 10 steady gait cycles, were used to minimize the objective function in (15). For the single-speed modes, the 10 steady gait cycles were all from one speed, while for the inter-speed mode, the 10 gait cycles were composed of two steady cycles from each of five speeds. Similarly, by manually setting the muscle activation allocation gain δ i to be 1 and 0, we could also conduct the calibration procedures of the sEMG-driven and US imaging-driven HNMs. After the HNM calibration with different modes, new data sets (five gait cycles not involved in the calibration procedures) from walking trials at different speeds were used as input of the calibrated HNMs to predict net PF moment. The prediction results were evaluated through the comparison to the measured net PF moment from ID.

Statistical analysis
The validation procedures comprised three tests to assess the sEMG-US imaging-driven, sEMG-driven, and US imaging-driven HNMs' calibration and prediction abilities of the net ankle joint PF moment during the stance phase. In each type of HNM test, the root mean square error (RMSE), RMSE normalized to individual peak net PF moment ( N − RMSE ), RMSE normalized to individual body mass ( BM − RMSE ), and coefficient of determination ( R 2 ) values between the calibrated/predicted net PF moment and ID-calculated net PF moment were calculated to evaluate the calibration/prediction performance. The rational of using these three tests is given here. From a mathematical perspective, the proposed sEMG-US imaging-driven HNM is a model-based regression problem. Given that the corresponding data were all collected continuously and synchronously, and the ground truth from the ID calculation and the prediction from different HNMs are all continuous, the most intuitive metric to evaluate the prediction accuracy is to use the RMSE between the ground truth and prediction. However, due to the weight, height, and walking pattern's variations among subjects, the capability (peak value) of net plantarflexion moment generation during treadmill walking at the same speed is different, which would affect the prediction RMSE directly. Therefore, the direct RMSE values comparison among subjects is likely to introduce high deviations. An effective way to reduce the subjective variation of the prediction RMSE is to take the normalization to a subjective characteristic, like the body mass and peak net plantarflexion moment as we selected. The relationship between N − RMSE and BM − RMSE is that in the normalization calculation, they have same numerators but different denominators for the same subject. Although N − RMSE and BM − RMSE might provide potential redundancy information as can be seen from subsequent calibration and prediction summary results, they show different sensitivities to the walking speeds ( BM − RMSE is more sensitive). Another useful metric to evaluate the prediction accuracy is the coefficient of determination, known as R 2 , which ranges from 0 to 1 and measures the proportion of variation in the data that is accounted for in the model. In the current study, R 2 is used as a parallel evaluation metric as N-RMSE and BM-RMSE, and it is relatively not sensitive to the walking speeds. Although R 2 and N − RMSE might provide potential redundancy information, it is not always true that lower N − RMSE corresponds to higher R 2 , so we keep both to provide potential supplementary information.
Shapiro-Wilk parametric hypothesis test was used to determine the normality of the corresponding N − RMSE , BM − RMSE , and R 2 values of each calibration step under the single-speed modes and inter-speed mode, and each prediction step. According to the results of the Shapiro-Wilk test, two-way repeated-measure analysis of variance (ANOVA) or Friedman's tests followed by a Tukey's honestly significant difference tests (Tukey's HSD) was applied to evaluate the effect of HNMs' category and walking speed on different evaluation criteria, including N − RMSE , BM − RMSE , and R 2 values, during both calibration and prediction procedures. The significant difference level was chosen as p < 0.05 for all statistical tests. Effect sizes were reported as η 2 p and Cohen's d for main effects from ANOVA or Friedman's tests and pairwise comparisons from Tukey's HSD, respectively.

Ankle joint kinematics and kinetics, and plantarflexors neuromuscular features
The treadmill walking speed's effect on human ankle joint kinematics, kinematics, and lower extremities' muscles activities has been extensively investigated and discussed in a recent study [57], but without results related to the architectural change of those muscles. In this study, the time sequences of both LGS and SOL's muscle thickness during the recorded walking duration were extracted by using UltraTrack. Take the SOL muscle of Participant Sub08 as an example, Fig. 3 shows the MT tracking results throughout the 20 s walking duration at each walking speed. The red solid and blue dashed curves represent the MT change of the SOL muscle with and without key-frame correction, respectively. The results demonstrate that the key-frame correction could significantly reduce the time-related drift, which is due to the tracking error accumulation along with the walking duration (exampled US imaging video with the LGS and SOL's MT tracking results at multiple walking speeds can be referred to Additional files 2, 3, 4, 5, 6). The reported results here are consistent with the studies mentioned in [55].
The continuous results of the ankle joint kinematics and kinetics throughout the recorded 20 seconds at each speed were segmented as a percent of gait cycle from 0 to 100% according to the GRF measurements, which are shown in Fig. 4. The solid lines and light shadowed areas report the mean and one standard deviation (SD) values of the right ankle joint angular position, velocity, and net PF moment changes across all gait cycles at each walking speed on Sub01 (data from other participants can be found in the Additional file 1: Fig. S1-S5), where 0% and 100% represent the time instants when the heel-strike occurred in the current and consecutive gait cycles, respectively. By comparing the ankle joint moment calculated from ID, the averaged peak net PF moment across all walking gait cycles within each walking trial monotonically increased with the increase of walking speeds on each participant, indicating more power was needed for faster walking speed. The continuous results of neuromuscular features throughout the recorded walking duration were segmented as a percent of stance cycle from 0 to 100%, which are shown in Fig. 5. The solid lines and light shadowed areas report the mean and one SD values of the processed sEMG signals and US imagingderived MT signals from the LGS and SOL muscles on the right leg during multiple stance phase cycles at each walking speed on Sub01 (data from other participants can be found in the Additional file 1: Fig. S11-S15), where 0% and 100% represent the time instants when the heel-strike and toe-off occurred in the same gait cycle, respectively. The upper two rows present the LGS and SOL MT changes during the stance phase, while the lower two rows present the muscle's low-pass filtered sEMG changes. The results show small variations of MT and sEMG for both LGS and SOL muscles at heel-strike  across all walking speeds, but relatively high variations of sEMG for both LGS and SOL muscles at toe-off, where the sEMG signals increase with increased walking speed. For all five walking speeds, we observe that MT of LGS or SOL muscle is almost identical at heel-strike and toeoff, however, the processed sEMG signals of LGS or SOL muscle show higher value at toe-off than at heel-strike. By comparing the shadowed areas of the same feature across speeds, we observed that features' deviations are higher at a slower speed, like 0.50 m/s, than those at higher speeds. This implies that keeping a steady muscle contraction pattern at a slower speed is more difficult than at a higher speed.

Correlation analysis results
Overall, in Fig. 5, both temporal MT and sEMG signals from LGS and SOL muscles show a fairly strong correlation with the net ankle joint PF moment calculated from ID. For each participant, it was assumed the step-to-step variation is negligible, so we segmented all periods of the ankle joint net plantarflexion moment during the walking stance phase at each walking speed, and then calculated a correlation coefficient value between the net plantarflexion moment time sequence and each neuromuscular feature time sequence during the same stance cycle. Therefore, for a specific neuromuscular feature and walking speed, if the stance cycles are n 1 , n 2 , ..., n 10 for all ten participants, the number of correlation coefficients would be 10 i=1 n i . To address the person-to-person variation and to get each entry in Table 1, we applied the Fisher transformation to yield approximately normally distributed correlation coefficients. After the Fisher transformation, we calculated the mean and one SD values of the correlation coefficients for each walking speed and the neuromuscular feature, as listed in Table 1. The results show that a positive relationship exists between the processed LGS or SOL's sEMG signal and net PF moment, and between the LGS or SOL's MT and net PF moment. All the mean values of correlation coefficients are higher than 0.8 except for the SOL's sEMG signal at 1.50 m/s walking speed. From the results of a two-way repeated-measure ANOVA, we did not observe a significant difference in the correlation coefficient values among the four neuromuscular features ( p = 0.201 ) or among the walking speed ( p = 0.112 ), which indicates the MT and sEMG signals might have comparable capabilities to predict net PF moment across speeds.

Results of HNM calibration and net PF moment prediction
In the HNM calibration procedures shown in Fig. 2b, one essential input component is the muscle activation levels of both LGS and SOL muscles, either only based on sEMG signals, US imaging-derived MT signals, or data fusion between them. The representative demonstration of three types of muscle activation levels of both LGS and SOL muscles during the walking stance phase at 0.75 m/s on Participant Sub03 are shown in Fig. 6. The left and right column subplots represent the time sequence points of muscle activation levels that are used in the model calibration and prediction procedures. It is observed that the data fusion could balance the muscle activation levels from both sEMG signals and MT signals and further compensate for the activation level drift from only MT signals, which is potentially beneficial for the accuracy improvement of the net plantarflexion moment prediction. Fig. 7 presents the representative results of HNMs calibrations on Sub05 with data collected from each individual walking speed, where each solid centered line and the shadowed area represent the mean and one SD values of the measured and calibrated net ankle joint PF moment, and each curve is composed of 10 stance cycles at each walking speed. For this participant, the mean calibration RMSE values are 8.70 N · m, 9.10 N · m, and 7.90 N · m by using sEMG-, US imaging-, and sEMG-US imaging-driven HNMs under walking speed of 0.50 m/s, respectively. Results from other walking speeds also showed that the sEMG-US imaging-driven HNM could effectively reduce the calibration RMSE under either single-speed modes or inter-speed mode. To quantitatively evaluate the calibration performance by using different HNM categories and under different speed modes, the mean RMSE and R 2 values between the benchmark moment values from ID and calibrated moment values on each participant are summarized in Table 2. Based on the observation from Fig. 7, it appears that the relative error may be inflated during specific areas of the stance cycles. The current outcomes (i.e., N − RMSE , BM − RMSE , and R 2 values) may not be directly sensitive to those specific areas, because the three evaluation metrics report an averaged performance across the entire stance cycle. This limitation of the current outcome metrics can be addressed in future work by identifying systematic changes in the prediction accuracy throughout the stance phase cycle rather than across the entire cycle.
Overall, the calibration RMSE/R 2 values by using the sEMG-US imaging-driven HNM are lower/higher than those by using the sEMG-driven and US imaging-driven HNMs, respectively, for both single-speed modes and inter-speed mode on all participants. Additionally, 97.5 % of R 2 values are higher than 0.80, which indicates a strong linear correlation between the calibrated and ID-calculated net PF moment values. After calculating the N − RMSE and BM − RMSE values, the summarized calibration performance across all participants under both single-speed modes and inter-speed mode is presented in Fig. 8 Fig. 8, where letters and asterisks represent the statistically significant difference.
To verify the net PF moment prediction performance of the proposed sEMG-US imaging-driven HNM and compare the new HNM to both sEMG-and US imaging-driven HNMs, the three calibrated HNM categories, shown in Fig. 8, were applied to walking scenarios of five different speeds. Similarly, BM − RMSE , N − RMSE , and R 2 values between the predicted and ID-calculated net PF moments were evaluated and compared. Due to the space limitation, the detailed RMSE and R 2 values at each walking speed on individual participant in prediction can be found in Additional file 1: Tables S2 to Table S4. The summarized results in Fig. 9 show the mean and one SD of BM − RMSE , N − RMSE , and R 2 values in prediction at each of five speeds by applying three HNM categories. The subplots from the first column to the fifth column represent the applied calibrated HNMs were from five single-speed modes, while the subplots from the sixth column represent the inter-speed mode. In each subplot, the x−axis labels represent the prediction scenarios at five different speeds. Table 3 summaries the mean values of BM − RMSE , N − RMSE , and R 2 in the prediction by using the proposed sEMG-US imaging-driven HNM, where the bold numbers represent the cases that the calibration speed mode is the same as the prediction speed scenario, which always results in the lowest RMSE and highest R 2 values. From Fig. 9 and Table 3

Discussions
For the first time, this study investigated net ankle joint PF moment continuous prediction during the walking stance phase based on an HNM approach that combines both sEMG-induced and US imaging-induced muscle activation components. We defined this approach as sEMG-US imaging-driven HNM for the net PF moment prediction. The modified HNM was validated with quasiperiodic data, including sEMG signals and US imagingderived MT from plantarflexors, ankle joint kinematics, and kinetics measurements from multiple treadmill walking stance phases gait cycles at five different speeds. We focused on two plantarflexor muscles, e.g., LGS and SOL, when characterizing the PF function during walking. The basic objective was to establish a modified HNM for either LGS or SOL muscles, including (1) a weighted muscle activation model component, (2) a muscle-tendon unit geometry model component, and (3) a muscle dynamic contraction model component. This study also Fig. 7 Calibration results by using three HNM categories on Participant Sub05 under single-speed modes. The calibration performance by using sEMG-, US imaging-, and sEMG-US imaging-driven HNMs under 5 single-speed modes on Sub05. The red lines and light shadowed areas represent the mean and SD of the net PF moment from the ID calculation, the blue, green, black lines and light shadowed areas represent the mean and SD of the net PF moment from calibrations by using sEMG-driven HNM, US imaging-driven HNM, and sEMG-US imaging-driven HNM, respectively. a-e represent the results at each walking speed from 0.50 to 1.50 m/s    Fig. 8 Calibration performance across all participants by using three HNM categories. The mean and SD of BM − RMSE , N − RMSE , and R 2 values across all participants under both single-speed modes and inter-speed mode are displayed in upper, middle, and lower plots, respectively, during calibration procedures. Blue, red, and yellow bars represent the results by using sEMG-, US imaging-, and sEMG-US imaging-driven HNMs, respectively. Letters and asterisks represent the statistically significant difference from the two-way repeated-measure ANOVA followed by a Tukey's HSD, where letters stand for the evaluation of the walking speed factor (significant difference only exists between each two bar groups out of the left six groups without any overlap of letters) and asterisks stand for the evaluation of the applied HNM type factor (significant difference exists between each two bars out of the most right group when the asterisk is there) investigated the effects of HNM categories (sEMG-, US imaging-, and sEMG-US imaging-driven) and calibration speed modes (single-speed modes and inter-speed mode) on the net PF moment prediction performance during the walking stance phase. The ID-calculated net PF moment worked as the benchmark in both calibration and prediction. The results support our hypothesis that compared to sEMG-driven and US imaging-driven HNMs, the sEMG-US imaging-driven HNM accurately predicts the net PF moment. Further, the HNMs, when calibrated with the inter-speed mode, robustly predict the net PF moment. The evaluation of the proposed HNMs' calibration and prediction performance was based on three criteria: BM − RMSE , N − RMSE , and R 2 . In general, as reported in [58], the results were considered excellent if the N − RMSE values were smaller than 15%. For the calibration, the results in Fig. 7 presented a good calibration performance by using each of three HNMs on Sub05. Furthermore, from results across participants shown in Fig. 8, we observed that the mean N − RMSE values were all less than 10.92% regardless of the calibration speed modes and applied HNM categories. Primarily they were less than 8.30% irrespective of the calibration speed modes when applying the proposed sEMG-US imagingdriven HNM. For the prediction, as shown in Fig. 9 and Table 3, by applying the inter-speed mode calibrated sEMG-US imaging-driven HNM, the N − RMSE mean values were all less than 11.40% throughout those five walking speeds scenarios, which validated the prediction results were all excellent. Therefore, both calibration and prediction results present consistent findings that the N − RMSE values can be constrained within a satisfactorily small region by using the proposed sEMG-US imaging-driven HNM with the inter-speed mode calibration.
The superior net PF torque prediction of the sEMG-US imaging-driven HNM can be mainly due to the reason that the fused signals capture complementary mechanical and neural aspects of the muscle contraction [25]. information [25]. Specifically, sEMG signals measure electric potentials generated by muscle motor units. The amplitude and density of sEMG signals linearly correlate with the number of firing neurons, which offers a physical measurement of the microphysiological response [59]. US imaging directly visualizes muscle change at the macrophysiological performance [60], when the same group of neural motor units is activated. Thus, sEMG features and US imaging features provide information from electrical and mechanical aspects, respectively, in response to the same physiological stimulus. As a preliminary study of the sEMG-US imagingdriven HNM for walking on a treadmill across multiple speeds, the results are promising and can help overcome the challenges of joint motion intent detection in volitional control of assistive devices. However, there are still some limitations in the current study. First of all, only a few US imaging-derived parameters, MT from both LGS and SOL, measured muscle activation levels indirectly (a macrophysiological perspective). Other US-derived parameters, like FL and PA [21,25,27,35], of the MTU geometry model may further enrich the neuromuscular model. A potential difficulty here is the inability to capture a larger region of interest that contains entire muscle fascicles, mainly due to the small dimension of our US transducer. Another difficulty is to continuously track the muscle fascicles of the plantarflexors within the visualized region of interest during dynamic tasks, like walking. If always visualized, more geometry parameters from US imaging signals will likely further improve the HNM accuracy. Second, the US imaging-derived muscle activation was defined as a normalization function in (4). We also assumed a linear monotonic relationship between the MT of LGS or SOL and net PF moment change during the stance phase in Table 1. However, in [29,61], it was observed that despite an increase in MT of the LGS and SOL muscles during the stance phase, their muscle activations decreased in the last 10 to 20% of the stance phase. Also, in fast walking, MT of SOL was observed to decrease in the first 10% of stance, while the muscle activation increased. The above inconsistency may result from the difference in experimental setup, sEMG sensor or US transducer placement on the targeted muscles, MT tracking approach from US imaging signals, etc. In future work, a further comparison will be conducted to evaluate the muscle geometry parameters change and sEMG signals change during the same walking tasks. As detailed in the pre-processing data section, a key-frame correction [55] was applied to reduce the gait cycle-related drift of MT tracking results significantly. However, a subjective error in the manual operation is inevitable. Therefore, it was not practical to altogether remove the gait cyclerelated drift of US imaging-derived MT. This work has extended the preliminary results of voluntary isometric ankle PF studies [34,35] to dynamic walking studies on non-disabled participants. Further experiments on patients with impaired plantarflexors are needed to verify the practicability of the proposed sEMG-US imagingdriven HNM for net ankle joint PF moment prediction. We also observed that the calibrated HNMs with the inter-speed mode increased the net PF moment prediction robustness across different speed scenarios. The development of inter-speed HNM modeling is crucial because it will potentially allow a consistently accurate prediction performance of the net PF moment across multiple walking speeds. Furthermore, it will enable predicting and analyzing the dynamics of a more extensive set of walking speeds than was possible with a singlespeed mode, resulting in a deeper understanding of the neuromuscular dynamics during the human walking stance phase. Finally, it may open up to the development of robust neuromuscular human-machine interfaces for the simultaneous and proportional control of wearable assistive devices such as powered orthoses and prostheses under different walking speeds.
In this preliminary study, although the combination of US imaging and sEMG data in the model improved the prediction of the ankle moment, it was unknown if the prediction was better only because in this combination the solution space was higher by adding parameters in the optimization problem. To address this ambiguity, we have conducted additional model calibration and prediction procedures for comparison, where the time sequence of the US imaging-derived muscle activation was replaced by a randomly generated signal between 0 and 1. In this way, we could reproduce the same solution space just like the proposed sEMG-US imaging-driven HNM, so we would be able to investigate if the prediction benefits still exist or not. We took a representative participant as an example, and present the calibration performance and prediction performance by using five different HNMs, as shown in the Additional file 1: Fig. S18. It was observed that when the random signal (RS)-derived muscle activation levels of LGS and SOL muscles were introduced, although the solution space was increased, just like the proposed approach to fuse sEMG and US imaging signals, the calibration and prediction performance were both deteriorated when compared to the case where only sEMG-derived muscle activation levels were used. Therefore, based on the additional analyses, we believe that the better net plantarflexion moment prediction by using the proposed sEMG-US imaging-driven HNM was not only because in this combination the solution space was higher by adding parameters in the optimization problem, but more importantly, the US imaging-derived muscle thickness feature introduced complementary mechanical deformation characteristics during the muscle contraction, which dominantly resulted in the better prediction performance.

Limitation and future work
The key point to implement the proposed sEMG-US imaging-driven HNM is to find the allocation gain δ i for both LGS and SOL muscles through solving the nonlinear least square optimization problem. However, one possible limitation of this optimization process is that the allocation gains are speed-dependent and subjectdependent, which means that the optimal values of those two allocation gains vary along with the model calibrations when using data collected from different walking speeds on the same participant or from different participants at the same walking speed. So it is relatively challenging to generate a constant allocation gain for either LGS or SOL muscles across multiple speeds or general participants. The US imaging-derived muscle thickness signal is potentially sensitive to the position and orientation of the US transducer relative to the shank segment, which may affect the prediction accuracy of the proposed approach. In the current study, we did not track the position and orientation of the US transducer relative to the shank segment due to the assumption that the relative motion of the US transducer is negligible because a customized US transducer holder with an arc structure and elastic velcro straps were used to stabilize the position and orientation of the US transducer during the dynamic walking tasks. However, one possible limitation of this study is that the assumption may not hold for a faster walking speed. Therefore, a further study would be to investigate the effect of the relative position and orientation change between the US transducer and the shank segment to the prediction accuracy of the proposed sEMG-US imaging-driven HNM approach.
Another limitation of the current study is that we only considered the agonist LGS and SOL muscles to predict the ankle joint net PF moment, and the effects of the MGS and TA muscles were not considered. Although in the experiments we collected the sEMG signals from LGS, MGS, SOL, and TA muscles, ultrasound imaging data was collected from the LGS and SOL muscles only. Therefore, the fusion results were only based on these two muscles. The main reasons that we did not consider the effects from the MGS and TA muscles are (1) existing evidence suggests that sEMG signals from LGS and MGS are comparable [62], and similar results in the current study also can be seen in the Additional file 1: Fig.  S19; (2) during the walking stance phase, plantarflexor muscles (LGS and SOL) are the dominant muscles that generate the net ankle joint moment. Results of the processed sEMG signals showed that the activation of the TA muscle are relatively small during the walking stance phase compared to the plantarflexor muscles, as can be observed in the Additional file 1: Fig. S19. However, even though the TA muscle activation is low during the walking stance phase, the effect from the passive contraction from this antagonist muscle needs to be further investigated in future work.

Conclusion
This paper investigated the feasibility of an sEMG-US imaging-driven HNM to predict net ankle joint PF moment during the stance phase across multiple walking speeds. The results showed that on average, the net PF moment prediction RMSE, normalized to peak net PF moment/body mass, was significantly reduced when using the sEMG-US imaging-driven HNM, compared to sEMG-driven and US imaging-driven HNMs. The results also showed that the calibrated HNMs with the inter-speed mode were more robust for net PF moment prediction at different speed scenarios than HNMs with single-speed modes. The improved net ankle joint PF moment prediction during the dynamic walking tasks at different speeds could potentially lead to improved volitional control of assistive devices with more advanced and intelligent algorithms.
Additional file 1: Table S1. Participants' anthropometric characteristics and peak net PF moment during treadmill walking at each speed.  Figure S1. Right ankle joint kinematics and kinetics results from Participant Sub01. Figure S2. Right ankle joint kinematics and kinetics results from Participant Sub02. Figure S3. Right ankle joint kinematics and kinetics results from Participant Sub03. Figure S4. Right ankle joint kinematics and kinetics results from Participant Sub04. Figure S5.
Right ankle joint kinematics and kinetics results from Participant Sub05. Figure S6. sEMG raw data and processed data from both LGS and SOL muscles during 20 seconds walking experiments on Participant Sub01. Figure S7. sEMG raw data and processed data from both LGS and SOL muscles during 20 seconds walking experiments on Participant Sub02. Figure S8. sEMG raw data and processed data from both LGS and SOL muscles during 20 seconds walking experiments on Participant Sub03. Figure S9. sEMG raw data and processed data from both LGS and SOL muscles during 20 seconds walking experiments on Participant Sub04. Figure S10. sEMG raw data and processed data from both LGS and SOL muscles during 20 seconds walking experiments on Participant Sub05. Figure S11. US imaging-derived MT and processed sEMG signals during the walking stance phase from Sub01. Figure S12. US imaging-derived MT and processed sEMG signals during the walking stance phase from Sub02. Figure S13. US imaging-derived MT and processed sEMG signals during the walking stance phase from Sub03. Figure S14. US imagingderived MT and processed sEMG signals during the walking stance phase from Sub04. Figure S15. US imaging-derived MT and processed sEMG signals during the walking stance phase from Sub05. Figure S16.
The Ultratrack toolbox interface with US imaging-derived MT tracking procedures and results without key-frame correction on Participant Sub07 walking at 1.00 m/s. Figure S17.The Ultratrack toolbox interface with US imaging-derived MT tracking procedures and results with key-frame correction on Participant Sub07 walking at 1.00 m/s. Figure S18.


Effects of Nanobubbles in Dermal Delivery of Drugs and Cosmetics

Dermal delivery, which delivers drugs and cosmetics through the skin, has attracted significant attention due to its non-invasive and simple administration compared with oral or injectable administration. However, delivery of the ingredients through the skin barrier is difficult because the primary function of the skin is to protect the human body by preventing the invasion of contaminants. Although various techniques have been developed to overcome skin barriers, chemical toxicity, complicated processes, and expensive equipment still remain as obstacles. Moreover, green chemistry, which minimizes or eliminates the use of toxic chemicals, is required in the cosmetic industry. Thus, the development of a new method for dermal delivery is required. In this study, we provide a new method for dermal delivery using nanobubbles (NBs). NBs generated in oil improve the delivery effect of the active ingredients through the high Brownian motion and charge-balancing effect. Franz cell experiments and depigmentation experiments using the B16F10 melanoma cells were conducted to confirm the enhanced delivery effects. The system using NBs will contribute to the advancement of the dermal delivery of drugs and cosmetics.

Introduction
For thousands of years, people have applied various substances to the skin to achieve therapeutic effects [1]. Dermal delivery has received considerable attention in recent years in the fields of medicine and cosmetics because it can be a non-invasive alternative to conventional oral or injectable administration using the skin, which is the largest organ in the human body [1][2][3]. Indeed, the delivery of ingredients through the skin offers a variety of advantages, including the improved bioavailability of ingredients that suffer from the gastrointestinal environment or hepatic first-pass effects, the possibility of long-term drug delivery at a constant rate, reduced side effects, and improved patient compliance [4][5][6][7]. For cosmetics used in daily life to improve the skin condition, a non-invasive and self-administered method through the skin is preferred over a delivery method that is administered by professionals and accompanied by pain. However, the skin is the first physiological barrier that protects human and mammalian bodies against the invasion of pathogens and chemical/physical assaults [8,9]. This barrier function of the skin, which is a significant obstacle in dermal delivery, limits the delivery amount of the active ingredients [9,10]. Thus, the development of formulations (gels, hydrogels, organogels, nano emulsion, micelle, liposome, and exosome), or techniques using physical force (electrical voltage, ultrasound, micro-needles, thermal ablation, magnetophoresis, photomechanical waves, and electron beam irradiation) are being continuously sought to overcome these limitations to enhance the delivery effect [11][12][13][14][15]. Nevertheless, chemicals commonly used in formulations are biologically toxic and can cause skin irritation, and the equipment required for the technology using physical force is large and expensive [16]. Moreover, in the cosmetic industry, green chemistry is preferred or required, and the use of harmful chemicals should be minimized. To Nanomaterials 2022, 12, 3286 2 of 8 overcome these obstacles, new dermal delivery systems should be developed with minimal or no use of chemicals and equipment. Nanobubbles (NBs), defined as gas cavities with a submicron scale, have drawn increased attention in the field of drug delivery systems as a vehicle due to their several novel properties [17][18][19]. First, NBs are invisible to the naked eye because their diameter is smaller than the wavelength of visible light, and a NB solution is completely transparent [20,21]. Second, NBs are highly stable in solution for several months because of stronger Brownian motion than the buoyancy [22][23][24]. Third, NBs, which have a high surface area, can load the ingredients between the NBs and the solvent interface [25][26][27]. Finally, a burst of NBs via external stimuli such as sonication and physical force generates energy that enhances the delivery of the ingredients by increasing the diffusion effect [28][29][30][31][32]. NBs with these unique characteristics are used for the delivery of the active ingredients in various fields, and they have great advantages, especially in the cosmetic industry [33][34][35][36]. Since they are invisible in the solution, they can be applied to various cosmetics products, and their long-term stability ensures that the efficacy of functional cosmetics will be maintained over time. Although NBs have an attractive advantage in the field of dermal delivery and the cosmetics industry, their effects on improving the delivery of ingredients without specific stimulus have not been reported yet on the dermal delivery system.
Here, we report a new dermal delivery technique using NBs contained in oil to improve the delivery effect of active ingredients. The NB system is prepared by the simple mixing of the NBs and the active ingredients, and it does not require additional chemical or complicated experimental procedures. The high Brownian motion of the NBs is expected to increase the diffusion of the active ingredient, and the charge-balancing effect of NBs is expected to improve the movements of the active ingredient in the skin barrier. NBs were generated using a novel, previously patented device, and Nile red coloring was added to the NB oil for the additional experiments. [37,38] After the preparation of the solution, the delivery effect of the active ingredients was confirmed using the Franz cell method, which is proposed as a reliable alternative to animal testing for skin absorption testing [39]. Franz cells are divided into two compartments with an artificial membrane mimicking the skin, and the active ingredients were delivered through the membrane [40][41][42][43][44]. Finally, the improved effect of the active ingredients by increasing the delivery effect was confirmed through an in vitro experiment. Delivery test was further conducted using α-bisabolol, which is a representative depigmenting ingredient [45][46][47]. Aliquots were collected from the delivery experiment and then were incubated with the B16F10 melanoma cells to confirm the depigment effect.

Materials and Methods
The Material: Medium chain triglyceride (MCT) was supplied from New&New Co., LTD. (Chungnam, Korea). Oxygen was purchased from Daehan Chemical Co., Ltd. Nile red was purchased from the TCI (Tokyo, Japan). Mixed Cellulose Ester (MCE) membrane was purchased from Scilab Korea Co., LTD. (Seoul, Korea). Nanoparticle tracking analysis (NTA; NanoSight NS300, Malvern Instruments, Worcestershire, UK,). The Franz cell, which was produced on a custom order from Dail Scientific Trade, Inc. (Seoul, Korea) was used to analyze the delivery effect of active ingredients. Microplate reader (Tecan Trading AG, Zurich, Switzerland,) was used to measure to Nile red and melanin concentration. B16F10 melanoma cells were obtained from the Korean Cell Line Bank (Seoul, Korea). All cell reagents for in vitro studies such as phosphate-buffered saline solution, Dulbecco's modified Eagle's medium (DMEM), fetal bovine serum (FBS), and penicillin-streptomycin were all purchased from Sigma-Aldrich (St. Louis, MO, USA).
Generation of NBs in oil: NBs were produced by placing 1 L of MCT oil in a tank and continuously passing it through a custom-made NB generator using a diaphragm pump. During the first 5 min, pure air was injected, and the device was operated to degas other gases. The operation of the generator was stopped when the concentration of NBs reached 2.99 billion/mL. The number of NBs was measured using the NTA device. For consistent measurements, the same condition of laser wavelength, slider shutter, slider gain, and detect threshold was used as 532 nm, 1500, 330 and 5 respectively.
Franz cell experiments: MCE membrane was fitted between the donor and receptor chamber of the Franz cell. Receptor chamber was filled with the ethanol aqueous solution (40% v/v). Nile red solution (500 µM, 1 mL) was loaded to the donor chamber. Aliquots (100 µL) were extracted from the receptor chamber and the same volume of ethanol aqueous solution (40% v/v) was added after extraction. Collected aliquots were analyzed using microplate reader.
Biocompatibility analysis of NB oil: Cell studies were performed to evaluate the biocompatibility of the NB Oil using the B16F10 (melanoma cells) cell line. The cells with a density of 3 × 10 3 cells/well were cultured overnight in DMEM containing 10% FBS and 1% penicillin/streptomycin at 37 • C in a 5% CO 2 incubator. Thereafter, 100 µL NB oil was added to each well, which was filled with 900 µL DMEM containing 10% FBS, 1% penicillin/streptomycin, and 1% ethanol. After 24 h, the MTT assay was conducted using a microplate reader at 570 nm absorption.
Live and dead cell assay: B16F10 (melanoma cells) cell line with a density of 5 × 10 3 cells/well was cultured for 24 h in DMEM containing 10% FBS and 1% penicillin/streptomycin at 37 • C in a 5% CO 2 incubator. Thereafter, NB oil was added to each cell. After 24 h, the live/dead cell assay was conducted using confocal microscopy. The distribution of live/dead cells was visualized by using fluorescein diacetate (FDA) and propidium iodide (PI) staining methods, respectively.
Depigmentation test: The B16F10 cell line with a density of 1 × 10 3 cells/well were cultured overnight in DMEM containing 10% FBS, 1% penicillin/streptomycin, at 37 • C in a 5% CO 2 incubator. Thereafter, α-bisabolol or aliquots were added to the cell and incubated for 4 h. After refreshing the media, alpha-melanocyte stimulating hormone (α-MSH) with known concentration (1 µM) was added and incubated for 96 h. Over time, melanin concentration was measured using the microplate reader.
Depigmenting index calculation: The concentration of melanin measured with a microplate reader is converted to a depigmenting index using the formula below. Melanin (untreated) is the concentration of melanin obtained from B16F10 cells incubated without α-bisabolol, and melanin (treated) is the concentration of melanin obtained from B16F10 cells incubated with an aliquot containing α-bisabolol.

Results and Discussion
General bubbles disappear quickly in the solution because of the buoyancy, which is an upward force exerted by a fluid higher than the gravity force of the bubbles, which is a downward force. The low stability of the bubbles in the solution is a major challenge for the application of bubbles, which are non-toxic and can effectively deliver various active ingredients. Therefore, for the bubbles to retain in the liquid, Brownian motion must be greater than their buoyancy to counteract the upward force. Brownian motion is the random motion of particles suspended in a medium and is affected by several factors, including the viscosity of the medium, absolute temperature, and spherical particle size. The Stock-Einstein equation explains that the Brownian motion increases as the particle size decreases at the same temperature and medium. Therefore, the generation of nanometer-sized bubbles is important to retain the bubbles in a liquid for a long time.
In this study, a previously patented device was used to produce the NB oil. This device constitutes of oil tank, generator, oil pump and gas pump. As oil and air that were stored in the oil tank circulate through the device by the pump, the gas-liquid mixture continuously collides with the structure in the generator to produce the NBs. As shown in Figure 1a, the generated NB oil shows a different phenomenon compared to the pure oil under the green light laser irradiation. NBs dispersed in the oil scatter the light and Nanomaterials 2022, 12, 3286 4 of 8 lead to the Tyndall effect, whereas the pure oil does not show a special phenomenon. NBs were further confirmed by NTA to measure the size and concentration. The enerated NB oil was diluted for accurate analysis because the high density of the NBs leads to interference between the light scattered by the NBs during the measurements. The concentration and size of the NBs were 2.99 billion/mL and 145.9 nm (Figure 1b).
In this study, a previously patented device was used to produce the NB oil. This device constitutes of oil tank, generator, oil pump and gas pump. As oil and air that were stored in the oil tank circulate through the device by the pump, the gas-liquid mixture continuously collides with the structure in the generator to produce the NBs. As shown in Figure 1a, the generated NB oil shows a different phenomenon compared to the pure oil under the green light laser irradiation. NBs dispersed in the oil scatter the light and lead to the Tyndall effect, whereas the pure oil does not show a special phenomenon. NBs were further confirmed by NTA to measure the size and concentration. The enerated NB oil was diluted for accurate analysis because the high density of the NBs leads to interference between the light scattered by the NBs during the measurements. The concentration and size of the NBs were 2.99 billion/mL and 145.9 nm (Figure 1b). The use of animals in dermal delivery testing in the field of cosmetics has been forbidden in recent years, and various experiments are being replaced with non-animal methods [33]. In particular, the Franz cell experiments are proposed as a reliable alternative to animal testing for skin absorption testing [34][35][36][37][38]. In this experiment, the Franz cell was used to confirm the skin delivery effect of the active ingredients. The Franz cell consists of two compartments (donor and receptor chambers) and a membrane is fitted between the chambers to form a passageway for the active ingredient (Figure 2a). The active ingredient can be delivered to the receptor chamber only through the membrane. Nile red was chosen to observe the active ingredient delivery and the concentration can be measured by fluorescence. The NB oil containing the Nile red was loaded into the donor chamber, and the receptor chamber was filled with aqueous ethanol solution (40% v/v), which has a high solubility of Nile red, thus preventing the precipitation. Aliquots (100 μL) were extracted from the receptor chamber to confirm the amount of Nile red that was delivered, and that the same volume of the ethanol aqueous solution (40% v/v) was added to maintain the volume of solvent in the receptor chamber. The change in the concentration of Nile red according to the addition of the solvent was corrected through the calculation. Aliquots collected from 0.5, 1, 2, 4, 8, 12, and 24 h at room temperature were analyzed using a microplate reader, and those results were calculated to the concentration ( Figure  2b). Strikingly, a higher concentration of Nile red was observed in the presence of NBs than in the absence of NBs. After 24 h, the final concentration of Nile red was 2.06 mM in the absence of NBs and 3.00 mM in the presence of NBs, which was about 1.45 times higher (Table 1). This improved delivery effect is thought to be due to the high Brownian motion of the NBs, which increases the diffusion of the active ingredient and the charge balance effect of the NBs, which improves the movement of the active ingredient in the membrane. The use of animals in dermal delivery testing in the field of cosmetics has been forbidden in recent years, and various experiments are being replaced with non-animal methods [33]. In particular, the Franz cell experiments are proposed as a reliable alternative to animal testing for skin absorption testing [34][35][36][37][38]. In this experiment, the Franz cell was used to confirm the skin delivery effect of the active ingredients. The Franz cell consists of two compartments (donor and receptor chambers) and a membrane is fitted between the chambers to form a passageway for the active ingredient (Figure 2a). The active ingredient can be delivered to the receptor chamber only through the membrane. Nile red was chosen to observe the active ingredient delivery and the concentration can be measured by fluorescence. The NB oil containing the Nile red was loaded into the donor chamber, and the receptor chamber was filled with aqueous ethanol solution (40% v/v), which has a high solubility of Nile red, thus preventing the precipitation. Aliquots (100 µL) were extracted from the receptor chamber to confirm the amount of Nile red that was delivered, and that the same volume of the ethanol aqueous solution (40% v/v) was added to maintain the volume of solvent in the receptor chamber. The change in the concentration of Nile red according to the addition of the solvent was corrected through the calculation. Aliquots collected from 0.5, 1, 2, 4, 8, 12, and 24 h at room temperature were analyzed using a microplate reader, and those results were calculated to the concentration (Figure 2b). Strikingly, a higher concentration of Nile red was observed in the presence of NBs than in the absence of NBs. After 24 h, the final concentration of Nile red was 2.06 mM in the absence of NBs and 3.00 mM in the presence of NBs, which was about 1.45 times higher (Table 1). This improved delivery effect is thought to be due to the high Brownian motion of the NBs, which increases the diffusion of the active ingredient and the charge balance effect of the NBs, which improves the movement of the active ingredient in the membrane.  To further demonstrate the biocompatibility of our novel system, in vitro experiments of the pure NB oil and the NB oil containing α-bisabolol with the B16F10 melanoma cell were undertaken. A MTT assay was performed to analyze the cell viability after incubation for 24 h with the NB oil that contained the different concentrations of α-bisabolol. The developed method shows nontoxicity at all concentrations and in the pure oil ( Figure  3b). In addition, the B16F10 cells were incubated with the NB oil to observe the cell viability using confocal. After incubation, the cell line was treated with FDA (green) for staining live cells and PI (red) for staining dead cells. As shown in Figure 3a, it was confirmed that most of the cells were live in the presence of the NB oil, which was not significantly different from that in the absence of the NBs. This result suggests that the NB oil has a high biocompatibility due to the use of only the solvent and the active ingredients; this means that it can be a safe system for the dermal delivery of various ingredients. After successfully demonstrating the delivery experiments of Nile red, a delivery test was further conducted using α-bisabolol, which is a representative depigmenting ingredient, to inhibit melanin synthesis related to hyperpigmentation. The B16F10 melanoma cells were treated with aliquots that were collected from the receptor chamber after 12 h, and the melanin concentration was measured using the microplate readers to calculate the depigmenting index [48]. The B16F10 melanoma cells, one of the skin cancers, produce melanin when incubated with α-MSH. However, when α-bisabolol, which inhibits the enzymes involved in melanin production, is applied to the cells, melanin production is reduced. Therefore, the amount of delivered α-bisabolol is estimated by comparing the calculated depigmenting index. The depigmenting index with both aliquots collected in the presence of the NBs shows that the depigmenting index was significantly increased to 400% by the delivered α-bisabolol, whereas a relatively small increase in the depigmenting index to 120% was observed with the aliquots collected in the absence of the NBs (Figure 3c). The difference in this depigmenting index is because of the delivered amount of α-bisabolol, which suggests that the NBs enhance the delivery of the active ingredients without physical stimulation or complex processes.  To further demonstrate the biocompatibility of our novel system, in vitro experiments of the pure NB oil and the NB oil containing α-bisabolol with the B16F10 melanoma cell were undertaken. A MTT assay was performed to analyze the cell viability after incubation for 24 h with the NB oil that contained the different concentrations of α-bisabolol. The developed method shows nontoxicity at all concentrations and in the pure oil (Figure 3b). In addition, the B16F10 cells were incubated with the NB oil to observe the cell viability using confocal. After incubation, the cell line was treated with FDA (green) for staining live cells and PI (red) for staining dead cells. As shown in Figure 3a, it was confirmed that most of the cells were live in the presence of the NB oil, which was not significantly different from that in the absence of the NBs. This result suggests that the NB oil has a high biocompatibility due to the use of only the solvent and the active ingredients; this means that it can be a safe system for the dermal delivery of various ingredients. After successfully demonstrating the delivery experiments of Nile red, a delivery test was further conducted using α-bisabolol, which is a representative depigmenting ingredient, to inhibit melanin synthesis related to hyperpigmentation. The B16F10 melanoma cells were treated with aliquots that were collected from the receptor chamber after 12 h, and the melanin concentration was measured using the microplate readers to calculate the depigmenting index [48]. The B16F10 melanoma cells, one of the skin cancers, produce melanin when incubated with α-MSH. However, when α-bisabolol, which inhibits the enzymes involved in melanin production, is applied to the cells, melanin production is reduced. Therefore, the amount of delivered α-bisabolol is estimated by comparing the calculated depigmenting index. The depigmenting index with both aliquots collected in the presence of the NBs shows that the depigmenting index was significantly increased to 400% by the delivered α-bisabolol, whereas a relatively small increase in the depigmenting index to 120% was observed with the aliquots collected in the absence of the NBs (Figure 3c). The difference in this depigmenting index is because of the delivered amount of α-bisabolol, which suggests that the NBs enhance the delivery of the active ingredients without physical stimulation or complex processes.

Conclusions
Dermal delivery through the skin has been researched with significant intere medical and cosmetic fields because of its simple and non-invasive administrat successfully developed a method that can effectively deliver active ingredients complicated equipment or procedures using NB oil. The NBs with a concentration billion/mL and a size of 145.9 nm was used for experiments, and when the delive of the active ingredients was compared, it was confirmed that the amount of the d ingredients was higher in the presence of NBs than in the absence of NBs regar temperature. Moreover, it was confirmed that when an aliquot collected in the p of the NBs was added to the B16F10 melanoma cells, the depigmenting index in 3.3 times compare to when the aliquot collected in the absence of the NBs was These experimental results indicate that NBs have the potential as a novel techn the field of dermal delivery.

Conclusions
Dermal delivery through the skin has been researched with significant interest in the medical and cosmetic fields because of its simple and non-invasive administration. We successfully developed a method that can effectively deliver active ingredients without complicated equipment or procedures using NB oil. The NBs with a concentration of 2.99 billion/mL and a size of 145.9 nm was used for experiments, and when the delivery effect of the active ingredients was compared, it was confirmed that the amount of the delivered ingredients was higher in the presence of NBs than in the absence of NBs regardless of temperature. Moreover, it was confirmed that when an aliquot collected in the presence of the NBs was added to the B16F10 melanoma cells, the depigmenting index increased 3.3 times compare to when the aliquot collected in the absence of the NBs was added. These experimental results indicate that NBs have the potential as a novel technique in the field of dermal delivery.