Characterizing the limits of human stability during motion: perturbative experiment validates a model-based approach for the Sit-to-Stand task

TL;DR

This study introduces and experimentally validates the Stability Basin, a model-based method to predict fall risk during Sit-to-Stand movements, showing high accuracy in identifying stability limits and potential fall scenarios.

Contribution

The paper presents a novel computational approach to determine stability regions during motion and validates it through a perturbative experiment with human subjects.

Findings

SB predicts trial outcomes with over 90% accuracy for failures.

SB predicts successful trials with over 95% accuracy.

SB outperforms three other methods in estimating stable regions by over 45%.

Abstract

Falls affect a growing number of the population each year. Clinical methods to identify those at greatest risk for falls usually evaluate individuals while they perform specific motions such as balancing or Sit-to-Stand (STS). Unfortunately these techniques have been shown to have poor predictive power and are unable to identify the magnitude, direction, and timing of perturbations that can cause an individual to lose stability during motion. To address this limitation, the recently proposed Stability Basin (SB) aims to characterize the set of perturbations that will cause an individual to fall under a specific motor control strategy. The SB is defined as the set of configurations that do not lead to failure for an individual under their chosen control strategy. This paper presents a novel method to compute the SB and the first experimental validation of the SB with an 11-person perturbative STS experiment involving forwards or backwards pulls from a motor-driven cable. The individually-constructed SBs are used to identify when a trial fails, i.e., when an individual must switch control strategies (indicated by a step or sit) to recover from a perturbation. The constructed SBs correctly predict the outcome of trials where failure was observed with over 90% accuracy, and correctly predict the outcome of successful trials with over 95% accuracy. The SB was compared to three other methods and was found to estimate the stable region with over 45% more accuracy in all cases. This study demonstrates that SBs offer a novel model-based approach for quantifying stability during motion, which could be used in physical therapy for individuals at risk of falling.

Figures

Tables

Table 1. Table 1: Statistics collected from the perturbative STS experiment are summarized below. Note that the onset of failure (Step %STS and Sit %STS) often occurred after the defined trial end, which is discussed in Sec. 2.4.1 .

STS Strategy	Total	CP	Steps	Sits	Trial Time (s)	Pert. Onset %STS	Step %STS	Sit %STS
Natural	$326$	$197$	$42$	$29$	$1.25\pm 0.23$	$50.34\pm 9.94$	$107.27\pm 20.92$	$96.16\pm 15.35$
MT	$322$	$198$	$40$	$19$	$1.13\pm 0.21$	$46.58\pm 11.03$	$114.30\pm 32.66$	$99.38\pm 15.21$
QS	$300$	$196$	$39$	$29$	$2.23\pm 0.89$	$50.87\pm 16.71$	$95.04\pm 20.53$	$85.26\pm 20.14$

Table 2. Table 2: The accuracy of the SBs generated using the BFF+FB controller model are reported below. This table reports a tally of the correct predictions for each STS strategy and trial type across subjects.

STS Strategy	Successful Trials	Step Trials	Sit Trials
Natural	248/255 – 97.25%	41/42 – 97.62%	27/29 – 93.10%
MT	257/263 – 97.72%	33/40 – 82.50%	17/19 – 84.75%
QS	213/232 – 91.81%	36/39 – 92.31%	27/29 – 93.10%
Combined	718/750 – 95.73%	110/121 – 90.91%	71/77 – 92.21%

Table 3. Table 3: The accuracy of SBs generated using the BFF+FB controller model are reported for CP trials at three different force levels of the applied perturbation. The results are aggregated across subjects and the three tested STS control strategies.

CP Force Level	Successful Trials	Step Trials	Sit Trials
Low	170/176 – 96.59%	15/15 – 100%	4/5 – 80.00%
Medium	125/132 – 94.70%	34/38 – 89.47%	24/27 – 88.89%
High	80/85 – 94.12%	61/68 – 89.71%	43/45 – 95.56%

Table 4. Table 4: The accuracy of the SBs formed using the BFF+FB controller model were compared to SBs formed using other methods. The results of the comparison are reported below, and are aggregated across subjects and STS strategies.

Basin Type	Successful Trials	Step/Sit Trials	False Successful Predictions	False Step/Sit Predictions
BFF+FB	718/750 – 95.73%	181/198 – 91.41%	17/735 – 2.31%	32/213 – 15.02%
LQR	71/750 – 9.47%	197/198 – 99.49%	1/72 – 1.39%	679/876 – 77.51%
FF+FB	351/750 – 46.80%	195/198 – 98.48%	3/354 – 0.85%	399/594 – 67.17%
naive	125/750 – 16.67%	195/198 – 98.48%	3/128 – 2.34%	625/820 – 76.21%

Equations

$$\tilde{x}(\tau )=\left[\begin{array}{c}{\tilde{r}}_x(\tau )\\ {\tilde{v}}_x(\tau )\\ {\tilde{r}}_y(\tau )\\ {\tilde{v}}_y(\tau )\end{array}\right].$$

$$\tilde{u}(\tau ,\tilde{x}(\tau ))=\left[\begin{array}{c}{\tilde{u}}_x(\tau ,\tilde{x}(\tau ))\\ {\tilde{u}}_y(\tau ,\tilde{x}(\tau ))\end{array}\right].$$

$$\dot{\tilde{x}}(\tau )=\left[\begin{array}{c}{\tilde{v}}_x(\tau )\\ \frac{1}{m}\left({\tilde{u}}_x(\tau ,\tilde{x}(\tau ))+{\tilde{d}}_x(\tau )\right)\\ {\tilde{v}}_y(\tau )\\ \frac{1}{m}\left({\tilde{u}}_y(\tau ,\tilde{x}(\tau ))+{\tilde{d}}_y(\tau )\right)-g\end{array}\right],$$

$$x(t)=\left[\begin{array}{c}{\tilde{r}}_x(T\cdot t)\\ T\cdot {\tilde{v}}_x(T\cdot t)\\ {\tilde{r}}_y(T\cdot t)\\ T\cdot {\tilde{v}}_y(T\cdot t)\end{array}\right].$$

$$u^{\text{lb}}\le u\le u^{\text{ub}}$$

$$u^{\text{lb}}(t,x(t))$$

$$u^{\text{ub}}(t,x(t))$$

$$\underset{{\text{ff}}^{\text{lb}}(t),{\text{ff}}^{\text{ub}}(t),K(t)}{\min }$$

$$+{\left(u^{\text{ub}}(t,x_i(t))-u_i(t)\right)}^2$$

$$u^{\text{lb}}(t,x_i(t))\le u_i(t)$$

$$u^{\text{ub}}(t,x_i(t))\ge u_i(t)$$

$$Z=\left\{y\in {\mathbb{R}}^n|y=c+G\beta ,-1\le \beta \le 1\right\}$$

$$\underset{{\beta }^{(t)}\in \mathbb{R},\beta \in {\mathbb{R}}^p}{\min }$$

$$G^{(t)}\beta =x(t)-c^{(t)},$$

$$|\beta |\le {\beta }^{(t)}$$

$$u(t,x(t))=u_{ol}(t)-K(x(t)-\bar{x}(t))$$

$$u(t,x(t))=\text{ff}(t)-K(t)x(t)$$

$$\underset{\text{ff}(t),K(t)}{\min }$$

$$RMSE(\hat{y},y)=\sqrt{\frac{{\sum }_{i=1}^N({\hat{y}}_i-y_i{)}^2}{N}}$$

Full text

Keywords: biomechanics, mathematical model, locomotion, fall risk, reachability, feedback control

1 Introduction

Falls are the leading cause of injury in people over 75 years old, resulting in reduced quality of life, increased healthcare costs, and accident-related death [1, 2, 3]. If at-risk individuals can be identified prior to injury, the likelihood of falling can be significantly reduced through intervention such as physical therapy [4, 5]. Fall risk generally results from instability arising from neuromuscular deficiencies or external perturbations. Clinical assessments for identifying individuals who would benefit from preventative care are currently limited to questionnaires [6, 7] and non-perturbative motor assessments [8, 9, 10]. Self-reported information often has low reliability [11, 12], and current clinical motor performance tests have low fall-prediction rates, especially for active older adults [13, 14, 15, 16]. For widespread use, stability assessments must combine predictive power with minimal experimental and computation time.

Since older adults are more likely to fall while in motion [17], several studies suggest that quantifying dynamic stability may help identify biomechanical deficiencies associated with an increased risk of falling [18, 19, 20, 21]. Thus, a number of model-based methods have been developed to assess stability during walking [21]. Among these, variability measures [22] and the maximum Lyapunov exponent [23] ranked highest overall in validity. However, these metrics only characterize a subject’s ability to recover from small perturbations. Currently, it is unclear whether the most useful estimates of stability are provided by measuring an individual’s ability to recover from small perturbations or by computing the largest possible perturbation that they can withstand without switching control strategies or falling [21]. Unfortunately, these methods cannot yet be fully compared since a verified technique for computing the maximum perturbation a subject can recover from does not exist. Additionally, the best-performing model-based methods are limited to periodic motion.

Stability during non-rhythmic motion is of interest, especially since difficulty with aperiodic tasks such as Sit-to-Stand (STS) are strongly correlated with falls in older adults and necessary for maintaining independence and quality of life [24]. The motion of a person’s center of mass (COM) during STS can be modelled as an inverted pendulum, which requires an appropriate amount of angular momentum at seatoff to successfully complete the task [25]. Drawing from this idea, metrics of stability for STS are generally based on an individual’s initial COM velocity or acceleration [26, 27, 28]. However, a stability metric that considers data only at the onset of movement disregards valuable information about the control strategy used by the individual. In fact, human control strategies for STS lie on a spectrum ranging from Quasi-Static (QS), in which little momentum is used and the body position is statically stable throughout the motion, to Momentum-Transfer (MT), which is statically unstable and extensively uses momentum to achieve standing [29].

To account for the dynamic differences in STS motions under distinct control strategies, Shia et. al introduced the Stability Basin (SB), computed using individualized pendulum models of STS with linear-quadratic regulator (LQR) controllers. For this pendulum model, the SB at a given time is equal to the set of pendulum states that are able to successfully achieve standing under the specified controller. Though this SB under an LQR controller successfully distinguishes between less and more stable STS strategies, its ability to identify exactly the set of states that can arrive at standing without switching strategies is unverified.

To test whether the SB can accurately estimate an individual’s stability for a particular task, perturbations must be introduced in a way that causes the states of the individualized dynamic model to exit the stable region. An accurate prediction of instability corresponds to an experimentally observed failure of the individual’s control strategy.

Here we use an $11$-person perturbative STS experiment, described in Sec. 2.1.1, to validate stability predictions generated by the SB method. Study participants attempted to complete the STS motion with their natural, QS, and MT control strategies while subjected to forwards or backwards cable pulls, applied to their approximate COM. We use a dynamic model, presented in Sec. 2.2, and subject-specific controllers to form individualized SBs, then test whether each subject’s SB can accurately predict a switch in control strategies in order to recover from the cable pull, described in Sec. 2.4. We consider common failure modes, such as a step or sit, to be observed changes in STS control strategy [30]. In contrast to prior work that assumed that humans use LQR control during STS while constructing the SB, we use a data-driven method to create bounds on the control input for a particular subject and control strategy, described in Sec. 2.3. Across all subjects and control strategies, the proposed method correctly predicts unsuccessful trials over $90\%$ of the time, and correctly identifies successful trials over $95\%$ of the time in a leave-one-out assessment. When compared to three other methods of assessing STS stability, the SB method out-performed each alternative. This perturbative validation experiment demonstrates that the SB method is accurate and reliable, showing promise for integration into current healthcare infrastructure.

2 Methods

The goal of this work is to characterize stability during STS. By stability, we mean the magnitude, direction, and timing of perturbations under which a subject can successfully stand up without stepping, sitting back down, or switching control strategy. We hypothesize that SBs can accurately represent stability.

The SB can be understood intuitively as follows. Given a dynamical model of STS, kinematic observations of a subject performing STS can be used to form trajectories of the model. The SB defines a boundary between “successful” and “unsuccessful” trajectories of the model, where “unsuccessful” implies the subject will have to take a step, sit back down, or switch control strategy. The SB is introduced pictorially in Fig. 1, with a formal definition and computational details given in Sec. 2.3. An accurate SB should predict when a subject’s STS motion will be successful or unsuccessful based on model trajectories alone. To test the accuracy of SBs, we performed a perturbative STS experiment, occasionally inducing subjects to step or sit back down, and then compared SB predictions to experimentally observed outcomes.

This section describes our framework for computing and testing the accuracy of SBs. First, we collect kinematic observations of a subject performing a perturbative STS experiment (Sec. 2.1). Individualized biomechanical models of STS are constructed for each subject (Sec. 2.2). Strategy-specific input bounds are used to compute the SBs (Sec. 2.3). We test whether the individual and strategy-specific SBs correctly predict when a subject steps or sits down in response to perturbation, and when they are successful (Sec. 2.4). Finally, we perform a comparison to other controller models and a naive method for estimating stability (Sec. 2.5).

2.1 Perturbative STS experiment

2.1.1 Experimental protocol

Subjects began in a seated position on a stool with their arms crossed. The height of the stool was adjusted so that the subject’s thighs were parallel to the ground. Subjects practised standing up from the stool and were asked to find a comfortable foot position, which was then demarcated with a line of tape.

Subjects were instructed to perform three different STS control strategies: their natural strategy, a MT strategy, and a QS strategy as done by Shia et al. [31]. Subjects watched a demonstration of each strategy, and then practised each strategy a minimum of 10 times prior to data collection. The following set of treatments were applied to each STS strategy:

• •

Nominal trials: Subjects stood five times from a comfortable foot position using the specified control strategy.

• •

Foot shift (FS) perturbations: Subjects varied their foot placement in two-inch increments in the anterior-posterior direction from their original position, which were demarcated by taped lines on the ground. Subjects stood up at each increment, moving their feet backwards until their heels left the ground, and then forward until a strategy shift was observed. We noticed that subjects at extreme anterior foot positions attempted to stand up by lunging forward into a squatting position, and excluded these trials from our dataset since they represented a shift in strategy from the subject’s nominal behavior.

• •

Cable pull (CP) perturbations: Subjects returned their feet to their original position. A cable system was attached to the subject’s waist and connected to two high-torque motors. These motor driven cables applied impulses to the subject as they rose. The cable pulls were restricted to the anterior-posterior direction, and were applied either forward or backwards with variable timing and force. Specifically, three peak force levels – low, medium, and high – were calibrated to each subject. The low force level was designed to rarely induce stepping or sitting during STS, while the high force level was designed to induce stepping or sitting approximately half of the time. Six trials were taken at each force level for each STS strategy, with three pulling forwards and three pulling backwards, in random order. For more details, see Appendix A.

2.1.2 Kinematic observations

A 10-camera PhaseSpace motion capture system collected kinematic observations of 36 markers at 480 Hz. C-Motion’s Visual3D biomechanics software was used to fit body segment models to each subject’s data [32]. MATLAB was used for all subsequent analyses [33]. A 6th-order Butterworth filter with a cut-off frequency of 2 Hz was used to filter joint position trajectories.

We estimate the motion of the subject’s center of mass (COM) from these kinematic data. First, we use a 3-segment model to track the motion of the subject’s shank, thigh, and head-arms-torso segments in the sagittal plane. The approximate COM of each segment is computed using anthropometric data [34], and combined to find the trajectory of the total body’s COM position throughout STS. Then, we obtain COM velocity and acceleration trajectories by numerically differentiating the COM position trajectories.

2.1.3 Classifying trials as successful/unsuccessful

Two clear instances of a subject switching control strategies during STS are taking a step and sitting back down. We refer to the trials in which these occurred as steps and sits, and collectively label them unsuccessful. At times, we will also refer to these trials as failures. We define the onset of stepping or sitting rigorously in Sec. 2.4.1. If no step or sit is detected, the trial is considered successful. Note, this scheme can classify trials in which subjects rocked far onto their heels or toes to maintain balance as successful, so long as no step or sit occurs.

2.2 Dynamic modelling

To describe our dynamic STS model and trajectory formation, we adopt the following mathematical notation. Let $A\times B$ denote the Cartesian product of sets A and B. Let $\mathbb{R}^{n}$ be the $n$-dimensional Euclidean space, and let $\tau\in[0,T]\subset\mathbb{R}$ denote a time drawn from a time interval. Without loss of generality, we assume each STS trial starts at time $\tau=0$.

2.2.1 Telescoping Inverted Pendulum Model

We use a Telescoping Inverted Pendulum Model (TIPM, Fig. 2) to model each subject’s COM motion in the sagittal plane during each STS trial [35]. This model is capable of describing the sizable displacements of a subject’s COM in the horizontal and vertical directions that occur during STS. The TIPM consists of a point mass of mass $m\in\mathbb{R}$ representing the subject’s COM. Individualized TIPMs are constructed for each subject by setting $m$ as the subject’s mass.

We define the origin as the initial position of the ball of the subject’s foot. We let the subscript $(\cdot)_{x}$ denote a quantity in the anterior-posterior direction, and the subscript $(\cdot)_{y}$ denote a quantity in the vertical direction. At each time $\tau\in[0,T]$, we denote the following quantities of the TIPM as:

• •

Positions: $\tilde{r}_{x}(\tau),\tilde{r}_{y}(\tau)\in\mathbb{R}$

• •

Velocities: $\tilde{v}_{x}(\tau),\tilde{v}_{y}(\tau)\in\mathbb{R}$

• •

Accelerations: $\tilde{a}_{x}(\tau),\tilde{a}_{y}(\tau)\in\mathbb{R}$

• •

Inputs: $\tilde{u}_{x}(\tau),\tilde{u}_{y}(\tau)\in\mathbb{R}$

• •

Cable-pull forces: $\tilde{d}_{x}(\tau),\tilde{d}_{y}(\tau)\in\mathbb{R}$

Let $\tilde{x}:[0,T]\to X\subset\mathbb{R}^{4}$ be a state trajectory, so that $\tilde{x}(\tau)\in X$ is the model’s state at time $\tau\in[0,T]$:

[TABLE]

Let $\tilde{u}:[0,T]\times X\to U\subset\mathbb{R}^{2}$ be an input trajectory, so that $\tilde{u}(\tau,\tilde{x}(\tau))$ is the input at that time and state:

[TABLE]

The time derivative of $\tilde{x}(\tau)$ at time $\tau$, denoted $\dot{\tilde{x}}(\tau)$, can be written as:

[TABLE]

where $g$ is the gravitational acceleration $9.81m/s^{2}$.

2.2.2 Nondimensionalizing time

The time that it takes to complete STS varies from trial to trial. Details on how the start and end of each STS trial are chosen are provided in Appendix B. To make accurate comparisons across trials, we normalize each TIPM state trajectory $\tilde{x}$ by the trial’s length, so that each trajectory occurs over the interval $[0,1]$, which we refer to as 0 to 100 %STS. We introduce a unitless time variable $t$, related to $\tau$ by $t=\frac{\tau}{T}$. This dimensionless variable allows us to introduce a normalized state trajectory $x:[0,1]\to X$, where

[TABLE]

We scale accelerations and forces by a factor of $T^{2}$ following similar logic.

For the rest of this paper, we only consider normalized trajectories, unless explicitly stated. Let $r_{x},r_{y},v_{x}$ and $v_{y}$ refer to normalized position and velocity trajectories, and let $a_{x},a_{y},u_{x},u_{y}$, $d_{x}$ and $d_{y}$ refer to normalized acceleration and force trajectories.

2.3 Computing Stability Basins

Three elements are required to form the SB for each of an individual’s STS control strategies:

The TIPM dynamics (3). 2. 2.

Bounds on the TIPM control input denoted $u$, a lower bound on the control input denoted $u^{\mathrm{\textnormal{lb}}}$, and an upper bound on the control input denoted $u^{\mathrm{\textnormal{ub}}}$. Note that

[TABLE]

where the inequality is taken element-wise. 3. 3.

A standing set, $X_{T}$, that encapsulates all model states observed to correspond to successful standing.

Given these elements, we define the SB as follows:

Definition 1

The SB $\subset[0,1]\times X$ is a subset of times (0 - 100%STS) and model states from which a trajectory of the model obeying the dynamics (3) and input constraints (5) will arrive at the standing set $X_{T}$.

In the rest of this subsection, we explain how we compute the input bounds $u^{\mathrm{\textnormal{lb}}}$ and $u^{\mathrm{\textnormal{ub}}}$ and the standing set $X_{T}$ from data. Then, we use these objects to generate the SBs.

2.3.1 Bounding control inputs using observations

To form the SB, we must specify the TIPM’s control input at each time and state. Modelling human motion requires this controller to be estimated from observations. In many previous studies, researchers have modelled human controllers as an open-loop “feedforward” component plus a closed-loop “feedback” component [36, 37, 31], finding a single controller that best models the data. However, since the form of the controller employed by humans is unknown, it is difficult to say with confidence that a single controller adequately captures human behavior. Therefore, at each instance in time, we compute the range of control inputs of a certain form that can be applied while explaining the data we observed, thereby modelling a range of inputs a human may use to recover from perturbation. This allows us to compute SBs with a set of hypothetical controllers, instead of just one. Our method of SB computation, detailed in Sec. 2.3.3, is flexible to this type of controller specification, and in particular, we verify that this approach is not too conservative in Sec. 3.

We present this approach more formally in the following definition:

Definition 2

We model the inputs $u$ to the TIPM as a feedforward component plus a linear feedback component, where the feedforward component is drawn from a bounded set. We refer to this controller model as a bounded feedforward plus feedback (BFF+FB) controller. Specifically, given time $t$ and state $x(t)$, we define the bounds in (5) as

[TABLE]

where $\mathrm{\textnormal{ff}}^{\mathrm{\textnormal{lb}}}(t)$ and $\mathrm{\textnormal{ff}}^{\mathrm{\textnormal{ub}}}(t)\in\mathbb{R}^{2}$ are lower and upper bounds on the feedforward input, and $K(t)\in\mathbb{R}^{2\times 4}$ is a matrix of linear feedback gains.

We generate the parameters of the BFF+FB controller $\mathrm{\textnormal{ff}}^{\mathrm{\textnormal{lb}}},\mathrm{\textnormal{ff}}^{\mathrm{\textnormal{ub}}}$, and $K$ from data. This is done separately for each STS strategy. First, we compute the inputs for each successful trial within that strategy via inverse dynamics using (3) [38, Chapter 5]. Let $x_{i}(t)$ and $u_{i}(t,x_{i}(t))$ denote the state and input at time $t$ for the $i^{th}$ successful trial computed via inverse dynamics, and let $S$ denote the set of successful trials within a given strategy.

We solve a constrained linear least squares problem to obtain the bounds $u^{\mathrm{\textnormal{lb}}}(t,x(t))$ and $u^{\mathrm{\textnormal{ub}}}(t,x(t))$:

[TABLE]

where the inequalities are understood element-wise. This is a quadratic program, which can be solved rapidly and efficiently. The bounds produced by this program are depicted in Fig. 3.

Since many successful trials within an STS strategy involved cable-pull perturbations, we make the following remark:

Remark 3

When computing $u_{i}(t)$ via inverse dynamics for successful CP trials, we discard the portion of the trajectory during which the perturbation was active. This is possible because the perturbations were time synchronized with the motion capture data.

2.3.2 Generating the standing set

The goal of the STS motion is to arrive at a standing configuration, which we represent as the standing set $X_{T}$. To define $X_{T}$, we must first introduce a geometric object called a zonotope. A zonotope $Z$ is a polytope in $\mathbb{R}^{n}$ that is closed under linear maps and Minkowski sums [39], and is parameterized by its center $c\in\mathbb{R}^{n}$ and generators ($g^{(1)},...g^{(p)}$) shown in Fig. 4. We write the generators as columns of the generator matrix $G\in\mathbb{R}^{n\times p}$. A zonotope describes the set of points that can be written as the center $c$ plus a linear combination of the generator matrix columns, where each element of the coefficient vector $\beta\in\mathbb{R}^{p}$ must be between $-1$ and $1$:

[TABLE]

where the inequalities are applied element-wise and $\mathbf{1}$ is a vector of ones of the appropriate size.

We define $X_{T}\subset X$ as a zonotope that encompasses the final states (i.e., the states at 100%STS) of all successful trials observed for a given STS strategy. Specifically, $X_{T}$ is a zonotope with 4 generators computed from the final states as proposed by Stursberg [40, Section 3], where each generator is expanded by 5% to avoid observed states lying on the edge of the set. The standing set generation process is depicted in Fig. 4.

2.3.3 Computing SB via reachability analysis

The reader may notice that the SB (Defn. 1) is defined as a set of states through time that satisfy some dynamics specified by an ordinary differential equation (ODE). Generally, numerical integration can be used to find solutions of an ODE starting from a single point. As we are interested in finding all possible solutions of the ODE, subject to constraints, we use a method called reachability analysis to compute the SBs. Rather than flowing a single point through a vector field, reachability analysis flows a set through a vector field. The set of all states that can be attained by the system, starting from some initial set, is referred to as the forwards reachable set. Conversely, all states that can end in some final set are referred to as the backwards reachable set.

To compute the SBs, we use an open-source reachability toolbox called CORA, which efficiently handles difficulties associated with representing set dynamics [41, 42]. CORA represents the SB as a zonotope at each of a finite collection of time steps, which are subintervals of the interval $[0,1]$ of length $\Delta t$. With a slight abuse of notation, we let $t$ act as an index, so that $Z^{(t)}$ is the zonotope representing the SB over the time step that contains $t$. Briefly, CORA works by linearizing the dynamics at each time step about the center of $Z^{(t)}$, and obtaining $Z^{(t+\Delta t)}$ by multiplying $Z^{(t)}$ by an overapproximation of the matrix exponential over that time step. It then expands $Z^{(t+\Delta t)}$ to account for the effects of inputs and linearization error. Each SB was formed using a time step $\Delta t$ of 0.005 (i.e., 0.5%STS), so that 200 zonotopes represent the SB over the interval $[0,1]$.

CORA allows the input to the system at each time step to be drawn from a set, allowing us to use the BFF+FB controller model and let the control input be defined as in (6) and (7). CORA therefore assumes the input at a given time and state can take any value within the input bounds to create the SB. Since the SB is the set of states through time that arrive at the final standing set $X_{T}$, it can be computed as a backwards reachable set. In practice, this means we treat $X_{T}$ as an initial set, and use CORA to flow the set backwards in time under the negative of the dynamics (3).

2.4 Evaluating Stability Basin accuracy

We hypothesize that the SBs can predict when a subject must switch their control strategy in response to perturbation to avoid falling. Previously in Sec. 2.1.3, we defined strategy switches as stepping or sitting back down. In this subsection, we explain how we determine the onset of stepping or sitting. Then, we use the SBs to predict whether or not a strategy switch will occur by checking if an STS trajectory exits the SB at any point. Finally, we detail the evaluation procedure we used to test the accuracy of the SBs’ predictions.

2.4.1 Determining onset of failure

To test the predictive power of the SBs, we first identify the onset of failure during an STS movement, $t_{f}\in[0,1]$, and show that the SB predicts failure prior to that time. For step trials, $t_{f}$ is defined as the time when the subject moves the toes of either foot more than 3 inches in the anterior-posterior direction from its starting point, as measured by motion capture markers placed on the subject’s foot. Although previous studies measured ground reaction forces to identify the onset of steps [30], the distance threshold was used here to utilize motion capture data. For sit failures, $t_{f}$ is defined as the time when both $v_{x}$ and $v_{y}$ become negative, i.e. the COM begins to move back towards the seat. If a subject both stepped and sat down, $t_{f}$ is defined as the earlier of the two.

Since we are using a procedure based on an average nominal trial for aligning and segmenting trials (detailed in Appendix B), it is possible for the onset of stepping or sitting to occur after the trial end time $t=1$. The end times chosen by the segmentation procedure do not necessarily imply that the subject has reached standing, but are just a means to achieve a consistent segmentation of the data. We consider a trial with a step or sit occurring later than the trial’s end to be unsuccessful, and set $t_{f}=1$ in this case.

2.4.2 Checking trajectories

Given a trajectory $x$ from an STS trial, we use a linear program to determine whether that trajectory lies inside the SB at each time step. Recall that each zonotope $Z^{(t)}$ representing the SB is parameterized by its center $c^{(t)}\in\mathbb{R}^{4}$ and generator matrix $G^{(t)}\in\mathbb{R}^{4\times p}$. At each time $t$, we represent the state $x(t)$ as a linear combination of the generator matrix’s columns, and use a linear program to find the representation that has the smallest maximum coefficient $\beta^{(t)}\in\mathbb{R}$:

[TABLE]

where the absolute value and inequalities are applied element-wise. If $\beta^{(t)}$ is less than or equal to 1, the point $x(t)$ is inside the SB at that time step, depicted in Fig. 5.

2.4.3 Evaluation procedure

The SBs are used to predict the success or failure of every trial within each STS strategy. The accuracy of each SB is then tested by comparing its predictions to the experimentally observed outcomes of each STS trial.

Given the state trajectory $x$ observed for an STS trial, the linear programs described in (10) are used to predict success/failure. If a state trajectory $x$ remains inside the SB at each time step, the trial is predicted to be successful. If $x$ exits the SB at any time step, the trial is predicted to fail. To ensure that the SBs’ predictions are made fairly, we make the following restrictions:

For CP trials, the prediction uses only the portion of the SB after the onset of perturbation. 2. 2.

For trials in which failure was observed, the prediction uses only the portion of the SB before failure occurred (i.e., $t\leq t_{f}$ as in Sec. 2.4.1).

We perform separate types of evaluation for successful trials versus unsuccessful trials. For successful trials, we employ a leave-one-out procedure to avoid training and testing on the same data. After forming the standing set, we leave one successful trial out of the BFF+FB computations in (8). We then construct an SB using this BFF+FB controller to test the trial that was left out. We repeat this procedure until all successful trials have been tested.

For unsuccessful trials, we form a single SB using a BFF+FB controller formed using all observed successful trials. We use this single SB to predict the outcome of each unsuccessful trial within a given strategy.

2.5 Comparisons

We compare SBs formed using the BFF+FB controller model (Defn. 2) to two other controller models, as well as a naive method for estimating stability. The same validation procedure described in Sec. 2.4.3 is performed for each.

2.5.1 LQR controller

Previously, Shia et al. proposed to model STS using an LQR controller about a nominal STS trajectory [31]. To build this controller, we first form a single average nominal trajectory $\bar{x}$ for each of a subject’s STS control strategies by taking the mean of the five nominal STS trajectories. We generate an open loop controller, $u_{ol}:[0,1]\to U\subset\mathbb{R}^{2}$, for each average nominal trajectory via inverse dynamics. Then, we use LQR to design a linear feedback controller about the average nominal trajectory. We use a quadratic cost function specified by the state weighting matrix $Q=I_{4x4}$ and input weighting matrix $R=1x10^{-4}I_{2x2}$ to generate the feedback matrix $K$ [43, Chapter 16], where the weighting matrices were found empirically to produce the best results. At a given time and state, the input from the LQR controller can then be written as:

[TABLE]

2.5.2 FF + FB controller

We also test the efficacy of a traditional feedforward plus feedback controller (FF+FB), mentioned in 2.3.1. This controller differs from the BFF+FB controller because the feedforward component is not drawn from a set. Let the input $u(t,x(t))$ at some time $t$ and state $x(t)$ be described by the state feedback matrix $K(t)\in\mathbb{R}^{1\times 4}$ and a feedforward component $\mathrm{\textnormal{ff}}(t)$:

[TABLE]

To generate $K(t)$ and $\mathrm{\textnormal{ff}}(t)$ from data, we use linear least squares to solve a program very similar to (8):

[TABLE]

As in (8), $S$ represents the set of observed successful trials within a STS strategy, and $x_{i}(t)$ and $u_{i}(t)$ correspond to the observed states and input of the $i^{\text{th}}$ trial at time $t$.

2.5.3 Naive method

Stability can be estimated from observed perturbed trials by simply drawing a volume around the state trajectories of the observed successful trials. This method does not utilize reachability or a controller model, and relies on state trajectories alone. Here, we apply this method by generating a 4-dimensional zonotope at each time step that encloses all of the observed successful states. This is the same procedure that we employ to generate the standing set $X_{T}$ (described in Sec. 2.3.3 and illustrated in Fig. 4) applied at each instance in time over the course of the motion. Since this estimate of stability is constructed using zonotopes, we follow the same leave-one-out validation procedure described in Sec. 2.4.3.

3 Results

We collected 948 trials from 11 participants (three female and eight male; ages 18-32; height 1.70 $\pm$ 0.12 m; body mass 65.4 $\pm$ 10.2 kg), including 213 unperturbed, 194 foot-shifted, and 591 cable pull trials. Each subject gave their informed written consent, and had no physical or balance disorders which could affect their ability to perform STS. The subjects performed STS trajectories with generally the same characteristics as the (different) subjects examined by Shia [31]. Even though subjects were allowed to determine the speed with which each STS was performed, each subject’s natural, MT, and QS STS strategies showed distinct kinematic characteristics, depicted in Fig. 6. In particular, the maximum and minimum horizontal accelerations of the COM for the MT ([-2.1, 2.2] m/s2) and QS ([-0.8, 0.6] m/s2) strategies were significantly different from each other as well as from the natural strategy ([-1.3, 1.5] m/s2), depicted in Fig. 7, with p-values $<0.001$ when comparing any two strategies. While the duration of natural and MT trials were similar (average times of 1.25s and 1.13s), QS trials lasted longer (average 2.23s), due in part to the low observed COM accelerations. The SB shape and volume are highly dependent on the kinematic and temporal characteristics corresponding to the strategy selected, as in previous work [31].

All foot-shift trials were considered successful. When using the MT strategy, subjects were able to stand with their feet further forward (max. $7.82\pm 2.75$ inches) than when using QS (max. $1.64\pm 1.21$ inches). Subjects attempting the extreme anterior and posterior foot positions were either unable to initiate the STS movement, or violated the parameters of acceptable trials stated in Sec. 2.1.1.

Cable pull perturbations during STS frequently induced a step or sit to avoid falling. Out of 591 cable pull trials, we observed failure in 198, or 33.5%. Specifically, we observed 42 steps and 29 sits for the natural strategy, 40 steps and 19 sits for MT, and 39 steps and 29 sits for QS. Fewer sits were observed for the MT strategy than in both the natural and QS strategies, likely due to the larger forward momentum at seatoff for the MT strategy.

SBs were generated on a laptop computer with a 2.7 GHz Intel Core i7 processor. SBs computed using the BFF+FB controller take $0.76\pm 0.014$ seconds to compute. To assess the accuracy of each SB, we compared the predictions to the experimental observation for each trial according to the procedure described in Sec. 2.4.3. The aggregate prediction rates of the SBs formed using the BFF+FB controller for each STS strategy across subjects are given in Tab. 2. Across STS strategies, successful trials (Fig. 8(a)) remain inside the basin 95.73% of the time. Of the failure trials, 90.91% of steps are predicted prior to onset (Fig. 8(b)), while 92.21% of sits are predicted prior to onset (Fig. 8(c)). There was no consistent effect of perturbation force on the predictive power of the method, detailed in Tab. 3.

Overall, the SB failure predictions are over 90% accurate; only 17 out of 198 unsuccessful trials are incorrectly predicted to be successful. To understand these false successes, we measured the maximum Euclidean distance of each unsuccessful trajectory to its nearest successful trajectory. The false successes are 34.24% closer to their nearest successful trajectories than are the true failure trajectories. Additionally, the onset of failure with respect to the onset of perturbation for false successes averaged 9.47% earlier than in the true failures. Note that it is plausible that individuals switched control strategies for some subset of the trials (e.g., from QS to natural), which could affect the accuracy of the strategy-specific SBs. Since we only characterize trials that are sits or steps as failures, it is possible that an individual switched strategies to successfully stand but would not have reached standing with their original strategy.

We compared the accuracy of SBs computed using the BFF+FB controller to SBs formed using other methods (Tab. 4. SBs formed using the LQR controller proposed in previous work [31] correctly predict 9.47% of successful trials. SBs formed using a traditional FF+FB controller correctly predict the outcome of 46.80% of successful trials. Finally, a naive method for estimating stability yields an accuracy of 16.67% for successful trials. Though each of the LQR, FF+FB, and naive methods yield nearly 100% prediction accuracy for trials where failure occurred, the low successful trial prediction rates mean that few trials remain within each SB, implying a severe underapproximation of the true stable region. We also give the false successful prediction and false failure prediction rates of each method in Tab. 4, and display these rates for each strategy in Fig. 9. Each of the LQR, FF+FB, and naive approaches for computing SBs has a high false failure prediction rate, further showing that each underestimates the size of the true stable region.

4 Discussion

In theory, stability can be directly characterized by perturbing a subject in an infinite number of ways and drawing a volume around the observed successful trials. Because this approach is practically infeasible, several model-based metrics have been developed for assessing dynamic stability. Many, such as Floquet Multipliers and Lyapunov Exponents, are limited to analyzing small perturbations during periodic motion. Only a few metrics assess the largest perturbations that can be withstood during motion [21], or the stability of aperiodic motion. In this area, Fujimoto and Chou [44] constructed estimates of the stable region for STS at a single time based on the extrapolated center of mass concept. We build on these methods by leveraging novel data-driven controller models, reachability analysis, and subject-specific dynamics to estimate stability over the course of the STS motion. By performing a perturbative experiment, we demonstrate that our SB-based approach accurately predicts stability across a diverse set of body morphologies, perturbations, and STS control strategies.

Our data-driven control model consists of a feedforward (“predictive”) and a linear feedback (“reactive”) component. Although there is evidence that humans employ feedforward and feedback terms for control [36], the exact form of these controllers remains unclear. Linear feedback applied about the COM has been sufficient to explain observed walking [45] and standing [46] behavior. However, we found the LQR controller proposed by Shia et al. [31] mischaracterizes the size and shape of the stable region because it is not trained on any of the observed perturbed trials. The classical FF+FB controller produced false failures, implying that this underestimates the true stable region of movement. Instead, our BFF+FB controller model draws from a set of feedforward inputs, plus a linear feedback term, that encapsulates all inputs from observed successful trials, resulting in improved performance over models that return a single input. This improvement in the accuracy of stability estimation is an important step towards clinical deployment.

In addition to increased accuracy, the computationally-tractable, individualized, and time-varying aspects of the SB approach make it well suited for widespread deployment and integration with existing clinical methods. The rapid and automated computation of SBs from kinematics alone minimizes the temporal and financial investment in performing this analysis, especially if coupled with computer vision techniques for pose estimation from digital cameras [47]. This can enable longitudinal studies examining how stability changes over time, and whether a smaller SB is correlated with fall risk, to be performed both accurately and easily. Although we make use of perturbative data for computing SBs here, one advantage of our proposed method is that it only relies on data from successful trials, and does not require a subject to be perturbed to failure. Thus, the SB method shows promise for characterizing the stability of subjects already at a high risk of falling or injury, such as elderly or rehabilitating subjects. Additionally, the TIPM and BFF+FB models are computed for each individual, enabling examination of a wide range of body morphologies and STS control strategies. Furthermore, stability throughout a task can be studied by measuring the size and location of cross-sections of the SB at particular times. By combining the SB approach with other clinical techniques, such as electromyography, clinicians can identify muscles associated with unstable portions of motion to target preventative and rehabilitative care. Such detailed analyses will undoubtedly provide deeper insight into the mechanisms underlying motor control.

Our approach can also be used in real time to inform controller design for wearable robotic devices because it provides time-varying boundaries on the region inside which an individual is stable. If a prosthetic device or exoskeleton detects that its user has exited the SB and become unstable, it can adapt its control scheme to aid in recovery. Furthermore, just as current tuning of controllers for prosthetic devices and exoskeletons focuses on minimizing metabolic cost, [48, 49], wearable robotic controllers can be optimized to increase stability by maximizing the volume of the SBs for the combined robot-human system. Thus, the high time-varying accuracy and low investment in computation demonstrate that the SB approach has the potential to increase access of personalized preventative and rehabilitative mobility care and inform the design of robotic systems for stable movement.

Ethics

The experimental protocol was approved by the University of Michigan Health Sciences and Behavioral Sciences Institutional Review Board, eResearch ID: HUM00020554.

Data, code and materials

Code and STS data may be found here: \hrefhttps://github.com/pdholmes/STS_SBhttps://github.com/pdholmes/STS_SB

Competing interests

We declare we have no competing interests.

Authors’ contributions

PDH conducted the experiment, analysed the data, wrote SB computation and validation code, and drafted the manuscript; SMD analysed the data, wrote SB computation code, and drafted the manuscript; XYF constructed perturbation equipment, conducted the experiment, and edited the manuscript; TYM conducted the experiment, supervised data analysis, and edited the manuscript; RV conceived and directed the study, and edited the manuscript.

Acknowledgements

We thank Art Kuo for helping to collect this dataset, and Victor Shia for his code and advice. We thank the 11 anonymous subjects for their participation in the study.

Funding

This work is supported by the National Science Foundation, under CAREER Award 1751093 and Graduate Research Fellowship Program Grant No. 1256260 DGE.

Appendix A Generating cable pull perturbations

Two high-torque motors were attached to cable-pulley systems and placed in front of and behind the test platform. The cables were attached to a belt around the subject’s waist, and the height of the pulleys were adjusted such that the cables were horizontal when the subject was standing. The pulleys were located approximately 3-4 feet in front of and behind the subjects. The cables were only attached for the CP condition, and detached for the nominal and FS conditions.

A custom-written LabVIEW [50] program was used to control the motor torques. During the experiment, a low torque was constantly commanded in both motors to keep the cables from going slack. The torque level was balanced between the anterior and posterior cables to minimize force bias felt by the subject. Force sensors on the cables were synchronized with motion capture data, so that the timing of the perturbation application in relation to kinematic data was known. Perturbations were manually activated by the toggling of a handheld switch. The experiment operator applied perturbations to the subject at variable times following seatoff to gather a range of recovery responses over the course of STS.

Perturbations were active over a 250 ms period, during which the applied force would ramp up to a specified peak, and then quickly ramp down to a baseline. The peak force to be applied was determined before the experiment began. Specifically, three peak force levels – low, medium, and high – were calibrated to each subject. The low force level was designed to rarely induce stepping or sitting during STS, while the high force level was designed to induce stepping or sitting approximately half of the time. Stepping and sitting are more formally defined in Sec. 2.4.1, and we report the number of steps and sits observed for each force level in Sec. 3. These force levels were roughly chosen based on the subject’s height and weight, and adjusted manually during a pre-experiment test session. Six trials were taken at each force level for each STS strategy, with three pulling forwards and three pulling backwards, in random order.

Appendix B Determining start and end of STS

To segment the continuous kinematic STS data into individual trials, an appropriate start and end time for each trial must be chosen. A consistent segmentation method is imperative for making comparisons across trials. Ideally, the start and end of each trial is chosen so that important kinematic features of STS are aligned.

We developed the following segmentation procedure, displayed in Fig. 10, which is designed to align the peak horizontal accelerations of STS trials while being robust to differences in trial duration. Portions of the procedure utilize a root-mean-square error (RMSE), defined as:

[TABLE]

where $\hat{y}=\{\hat{y}_{i}\}_{i=1}^{N},\hat{y_{i}}\in\mathbb{R}$ and $y=\{y_{i}\}_{i=1}^{N},y_{i}\in\mathbb{R}$ are two sequences defined by sampling two different trajectories at $N$ times.

The procedure generates an average nominal trial from the subject’s five nominal trials, and then uses it as a template to segment the rest of the trials. We segment the natural and MT strategies’ trials using the following procedure, but use a separate procedure for QS trials, detailed later on in this section.

Manually segment each trial, starting several seconds before the subject begins to stand and ending several seconds after they have reached standing. 2. 2.

To create an average nominal trial:

1. (a)

Take the manual segmentations of the subject’s five nominal trials, and align their COM position trajectories such that the times at which the peak horizontal COM velocities occur are coincident. 2. (b)

Average the subject’s five nominal trials together to form a single average nominal COM position trajectory. 3. (c)

Numerically differentiate the average nominal COM position trajectories to obtain velocity and acceleration trajectories. 4. (d)

Choose the start time ($\tau=0$) of the average nominal trajectory as when the horizontal COM acceleration first exceeds $20\%$ of its max. Choose the end time ($\tau=T_{\mathrm{\textnormal{nom}}}$) as when the vertical COM position first exceeds $99\%$ of its max. 5. (e)

The average nominal trajectory is now defined over $[0,T_{\mathrm{\textnormal{nom}}}]$. Follow the scaling laws in Sec. 2.2.2 so that it is defined over $[0,1]$. 3. 3.

To segment each trial:

1. (a)

Iterate over times $T\in[T_{\mathrm{\textnormal{min}}},T_{\mathrm{\textnormal{max}}}]$, where the interval is discretized into 100 samples. Let $T_{\mathrm{\textnormal{min}}}=0.75s$ and $T_{\mathrm{\textnormal{max}}}=1.5T_{\mathrm{\textnormal{nom}}}$. 2. (b)

Rescale the average nominal trajectory to be defined over $[0,T]$. 3. (c)

Align the scaled average nominal trial and the current trial so that the time at which the peak horizontal COM accelerations occur are coincident. 4. (d)

The trajectories overlap in a window of length $T$ seconds. Compute the RMSE between the horizontal COM acceleration trajectories over this window. 5. (e)

Find $T$ that minimizes this RMSE (to accurately segment CP trials despite the effects of perturbation, we only minimize the pre-perturbation RMSE for CP trials).

Since the QS strategy utilizes minimal forward momentum, it lacks a distinct peak in its horizontal COM acceleration trajectories when compared to the natural and MT strategies. This is shown in Fig. 6. Therefore, we developed a separate segmentation procedure for the QS strategy, where Steps 1 and 2 remain the same as before:

To segment each QS trial:

1. (a)

Iterate over times $T\in[T_{\mathrm{\textnormal{min}}},T_{\mathrm{\textnormal{max}}}]$, where the interval is discretized into 100 samples. Let $T_{\mathrm{\textnormal{min}}}=0.75s$ and $T_{\mathrm{\textnormal{max}}}=1.5T_{\mathrm{\textnormal{nom}}}$. 2. (b)

Rescale the average nominal trajectory to be defined over $[0,T]$. 3. (c)

Manually align the scaled average nominal trial and the current trial so that the time at which the peak vertical COM velocities occur are coincident. 4. (d)

The trajectories overlap in a window of length $T$ seconds. Manually choose $T$ that minimizes the RMSE between the scaled average nominal vertical COM velocity trajectory and the current trial’s COM velocity trajectory. (Prioritize the pre-perturbation difference for CP trials, to accurately segment CP trials despite the effects of perturbation.)