Coupling relations underlying the production of speech articulator movements and their invariance to speech rate

Abstract

Since the seminal works of Bernstein (The coordination and regulation of movements. Pergamon Press, Oxford, 1967) several authors have supported the idea that, to produce a goal-oriented movement in general, and a movement of the organs responsible for the production of speech sounds in particular, individuals activate a set of coupling relations that coordinate the behavior of the elements of the motor system involved in the production of the target movement or sound. In order to characterize the configurations of the coupling relations underlying speech production articulator movements, we introduce an original method based on recurrence analysis. The method is validated through the analysis of simulated dynamical systems adapted to reproduce the features of speech gesture kinematics and it is applied to the analysis of speech articulator movements recorded in five German speakers during the production of labial and coronal plosive and fricative consonants at variable speech rates. We were able to show that the underlying coupling relations change systematically between labial and coronal consonants, but are not affected by speech rate, despite the presence of qualitative changes observed in the trajectory of the jaw at fast speech rate.

Appendix: Simulations details

We validate methods for detecting coupling in simulation experiments where we can control the direction and the level of coupling, as well as amplitude or temporal variability. Throughout this study, we are only considering unidirectional coupling. We are interested in the generalized synchronization of two dynamical systems X and Y, i.e., we do not use coupling scenarios like $\dot{x}=F(x)+\mu (y-x)$, which would easily lead to complete synchronization in our setups.

1.1 Uninterrupted simulations of structurally different systems

In order to obtain results comparable to those obtained in previous works based on joint recurrence analysis (e.g. Romano et al., 2007), we investigate a Van der Pol oscillator X and a Rössler system. In a first scenario, we let X drive Y. The systems are described by the equations:

$$\begin{aligned}&\dot{x}=x_{2}, \nonumber \\&{{\dot{x}}_{2}}=0.1\left( 1-x_{1}^{2} \right) {{x}_{2}}-{{\omega }^{2}}{{x}_{1}}+\nu {{\sigma }_{{{x}_{2}}}}\xi , \end{aligned}$$

(8)

and

$$\begin{aligned} {{\dot{y}}_{1}}= & {} -{{y}_{2}}-{{y}_{3}},\nonumber \\ {{\dot{y}}_{2}}= & {} {{y}_{1}}+0.15{{y}_{2}}+\mu \frac{{{\sigma }_{{{y}_{2}}}}}{{{\sigma }_{{{x}_{1}}}}}{{x}_{1}},\nonumber \\ {{\dot{y}}_{3}}= & {} ({{y}_{1}}-10){{y}_{3}}+0.2. \end{aligned}$$

(9)

Here, $\xi$ is Gaussian white noise with a standard deviation of 1, which is added to the driving system. The driven system features a coupling term. The amount of noise is given by $\upnu$ and the coupling strength is given by $\mu$. $\sigma$ is the standard deviation of the respective components, evaluated for a typical realization of the trajectories without noise or coupling ($v=0$, $\mu =0)$, used for normalization. $\omega$ is a frequency parameter, and for $\omega =1$ both systems have approximately the same frequency, depending on the levels of noise and coupling. Since values around 1.0 (e.g., $\omega \in \left[ {0.7,1.3} \right] )$ quickly produce phase synchronization even for low levels of coupling, we choose X having half of Y ’s frequency by $\omega =0.5$ in the setup of long trajectories. X then completes half a cycle on average in the temporal interval $\left[ {0,2\pi } \right]$ while Y completes one. In this way, the influence of coupling is rather local and leads to generalized synchronization. The systems are numerically integrated using the forward Euler method and a constant step size of $dt=0.01$ in the interval $\left[ {0,200} \right]$. In this interval, X completes around 16 cycles, Y around twice as many. The length of the interval is chosen so that the number of oscillations of the slower X trajectory roughly matches the number of gestures found in the real data investigated in this study (Fig. 11).

Fig. 11

A noisy Van der Pol system X (Eq. 8) driving a Rössler system Y (Eq. 9), $t\in \left[ {0,200} \right]$, $v=0.3$, $\mu =0.3$

A ”warmup” interval of the same length is used and discarded afterwards, which guarantees that the systems are not in a transient state. The remaining 20,000 time steps are downsampled to $N=2000$ points for further analysis.

For investigating the methods’ independence of the special structure of driver and driven systems, we also consider the case of X driven by Y. The equations read

$$\begin{aligned} \dot{x}_1= & {} x_2 , \nonumber \\ {{\dot{x}}_{2}}= & {} 0.1\left( 1-x_{1}^{2} \right) {{x}_{2}}-{{\omega }^{2}}{{x}_{1}}+\mu \frac{{{\sigma }_{{{x}_{2}}}}}{{{\sigma }_{{{y}_{1}}}}}{{y}_{1}}, \end{aligned}$$

(10)

and

$$\begin{aligned} {{\dot{y}}_{1}}= & {} -{{y}_{2}}-{{y}_{3}},\nonumber \\ {{\dot{y}}_{2}}= & {} {{y}_{1}}+0.15{{y}_{2}}+\nu {{\sigma }_{{{y}_{2}}}}\xi ,\nonumber \\ {{\dot{y}}_{3}}= & {} ({{y}_{1}}-10){{y}_{3}}+0.2. \end{aligned}$$

(11)

As above, a stochastic component $\xi$ is added to the driver and a coupling term is added to the driven system.

1.2 Multiple-runs simulations

In order to simulate the conditions found in speech data, we also evaluate the equations several times during short time-intervals. This means that for a certain number of repetitions (here: 30), we generate short trajectories of an approximate length of $2\pi$ under different initial conditions, cf. Fig. 12. For this setup, we can choose $\omega =1$ (same frequency), because the systems are not run long enough to get phase locked and we still obtain generalized synchronization. The initial conditions are drawn randomly from a set of typical realizations of the long trajectories subject to the constraint that the difference in phase across repetitions and between the two systems is smaller than $\pi /2$, so a certain temporal alignment is guaranteed. The systems are numerically integrated in the interval $\left[ {0,6\pi } \right]$ (using $dt=0.01$ as above), where they complete approximately three cycles. The portion of time corresponding to the second cycle is then extracted (where the beginning of the cycle is set at the time point where $x_{1}$ is closer to 0) and each of the 30 trajectories obtained is downsampled to 80 points (Fig. 12).

Fig. 12

15 cycles of a Van der Pol system X (Eq. 8) driving a Rössler system Y (Eq. 9),  $v=0.3,\, \gamma =0.3,\, \mu =0.1$

1.3 Sources of variability in the simulations

We add synthetic dynamical noise and temporal variability to the systems (Lancia et al. 2014). We distinguish between internal and observational variability. Internal noise is induced to the driver before coupling, such that it also affects the dynamics of the driven. Observational variability is applied to both systems afterwards and, therefore, does not interfere with the coupling relation, but will make it even harder to identify.

1.4 Internal dynamical noise

Internal dynamical noise is introduced to the driver by the additive white noise term $\xi$. The effect of the noise term is modulated by the parameter v. Due to the normalization, e.g. for $v=0.3$, the additive noise has a standard deviation 30% of the respective original signal’s standard deviation $\sigma$. Through the coupling term, intrinsic noise is propagated to the driven system.

1.5 Internal temporal variability

First we describe the general case of a continuous time-deformed signal, following Lucero (2005). The methods for the discretized signals are introduced subsequently. The signal $\tilde{x} \left( t \right)$, $t\in \left[ {0,T} \right]$, is obtained by the original signal $x\left( t \right)$, via a mapping function $H\left( t \right)$ introducing the required amount of temporal variability:

$$\begin{aligned} \tilde{x}\left( t \right) =x\left( {H\left( t \right) } \right) \end{aligned}$$

(12)

First, consider a mapping function which is the identity $H\left( t \right) =t$ or, equivalently, a signal $x\left( {H\left( t \right) } \right)$ moving at constant speed $dH\left( t \right) /dt=1$. Now, let the signal speed be perturbed by temporal noise

$$\begin{aligned} \frac{dH\left( t \right) }{dt}=1+\gamma \phi \left( t \right) , \end{aligned}$$

(13)

where $\phi \left( t \right)$ is a Gaussian process with zero mean, standard deviation of 1 and an autocorrelation function $exp\left( {-\left( {\theta t} \right) ^{2}} \right)$. The correlation length parameter $\theta$ is adapted such that it features approximately five extrema per cycle. $\gamma$ controls the amount of temporal variation. $\gamma \phi \left( t \right) >-1$ is required for a monotonic mapping $H\left( t \right)$. Otherwise, locally negative slopes occur, which means the new trajectory contains time-reversed portions of the original trajectory.

The mapping function $H\left( t \right)$ is obtained by integration:

$$\begin{aligned} H\left( t\right) =t+\gamma \mathop {\int }\nolimits _{0}^{t} \phi \left( s\right) ds. \end{aligned}$$

(14)

For a discrete time-series $x\left( n \right)$ with dimensionless time-scaling $\left( {n=0,1,\ldots ,N} \right)$, temporal variability is introduced as follows. A discrete Gaussian process $\phi \left( n \right)$, satisfying $\phi \left( 0 \right) =0$, is generated and the mapping function

$$\begin{aligned} H\left( n \right) =n+\gamma \mathop \sum \limits _{k=0}^n \phi \left( k \right) \end{aligned}$$

(15)

is computed, satisfying $H\left( 0 \right) =0$. In a second step, the mapping is normalized by

$$\begin{aligned} \tilde{H} \left( n \right) =H\left( n \right) \frac{N}{H\left( N \right) } \end{aligned}$$

(16)

which leads to $\tilde{H} \left( n \right) =N$. Since a mapping is generally not integer-valued ($\tilde{H} \left( n \right) \notin N)$, an interpolation of $x\left( 0 \right) ,x\left( 1 \right) ,\ldots ,x\left( N \right)$ has to be used to compute the values of the time-deformed signal

$$\begin{aligned} \tilde{x} \left( n \right) =x\left( {\tilde{H} \left( n \right) } \right) \end{aligned}$$

(17)

between the original trajectory’s discretization points $0,1,\ldots ,N$. For the multiple-runs simulations, the temporal deformation is applied to each run separately. The effect of temporal deformation is illustrated in Fig. 13.

Fig. 13

a Temporal deformation of the Van der Pol system (Eqs. 8, 10), dashed: original $x \left( n \right)$, solid: deformed $\tilde{x} \left( n \right) ,\gamma =0.3$. b Corresponding time-shift $\tilde{H} \left( n \right) -n$. c Quantification of temporal deformation: amplitude root-mean-square error $RMSE/\sigma$ vs. different values of $\gamma$ in Eq. (14) for 100 Van der Pol simulations each. E.g. $,\gamma =0.3$ introduces a mean amplitude error of $0.87\sigma$

To introduce internal temporal variability, we do not let $Y\left( \cdot \right)$ be driven by $X\left( \cdot \right)$, but by a temporally deformed signal $X\left( {H\left( \cdot \right) } \right)$, where the amount of variability is given by $\gamma$. The analysis of coupling relations is then performed on $X\left( {H\left( \cdot \right) } \right)$ and $Y\left( \cdot \right)$.

1.6 Observational temporal variability

We let X drive Y without any temporal deformation. Afterwards, both signals are deformed by different mappings, $H\left( \cdot \right)$ and $G\left( \cdot \right)$. The amount of temporal variability is given by $\gamma _x$ and $\gamma _y$, respectively. In this way we can also allow different amounts of noise in $X\left( {H\left( \cdot \right) } \right)$ and $Y\left( {G\left( \cdot \right) } \right)$. Then we compare the time-deformed signals $X\left( {H\left( \cdot \right) } \right)$ to $Y\left( {G\left( \cdot \right) } \right)$ for detecting the coupling relation.

1.7 Simulations parameters

Purely deterministic simulations: ($v=0$, $\gamma =0$, $\gamma _X =0$, $\gamma _Y =0)$; simulations with intrinsic noise and temporal variability: ($v=0.3$, $\gamma =0.3$,, $\gamma _X =0$, $\gamma _Y =0)$; uninterrupted simulations with intrinsic noise and variability and observational temporal variability: ($v=0.3$, $\gamma =0.3$, $\gamma _X =0.3$, $\gamma _Y =0.3)$.