1. Neuroscience
Download icon

Vocal and locomotor coordination develops in association with the autonomic nervous system

  1. Morgan L Gustison
  2. Jeremy I Borjon
  3. Daniel Y Takahashi  Is a corresponding author
  4. Asif A Ghazanfar  Is a corresponding author
  1. Princeton Neuroscience Institute, Princeton University, United States
  2. Princeton University, United States
Research Article
  • Cited 5
  • Views 809
  • Annotations
Cite this article as: eLife 2019;8:e41853 doi: 10.7554/eLife.41853

Abstract

In adult animals, movement and vocalizations are coordinated, sometimes facilitating, and at other times inhibiting, each other. What is missing is how these different domains of motor control become coordinated over the course of development. We investigated how postural-locomotor behaviors may influence vocal development, and the role played by physiological arousal during their interactions. Using infant marmoset monkeys, we densely sampled vocal, postural and locomotor behaviors and estimated arousal fluctuations from electrocardiographic measures of heart rate. We found that vocalizations matured sooner than postural and locomotor skills, and that vocal-locomotor coordination improved with age and during elevated arousal levels. These results suggest that postural-locomotor maturity is not required for vocal development to occur, and that infants gradually improve coordination between vocalizations and body movement through a process that may be facilitated by arousal level changes.

https://doi.org/10.7554/eLife.41853.001

Introduction

Vocal development is typically thought of as the adaptive coordination of the vocal apparatus (i.e., the lungs, larynx and the mouth) and the associated muscles and neural systems that influence its activity. In adult animals, however, vocal behavior does not occur in isolation of the operations of other motor systems. Studies investigating the real-time coordination between vocal and locomotor outputs show that some vocalizations can be produced concurrently with locomotor activity and/or postural changes, while others cannot (Suthers et al., 1972; Blumberg, 1992; Fusani et al., 1996; Williams, 2001; Wong and Waters, 2001; Holderied and von Helversen, 2003; Branchi et al., 2004; Cooper and Goller, 2004; Berg et al., 2013; Dalziell et al., 2013; Hoepfner and Goller, 2013; Ota et al., 2015; Alves et al., 2016; Laplagne and Elías Costa, 2016; Ullrich et al., 2016). For example, the ‘A’, ‘B’ and ‘D’ song types of the lyrebird (Menura novaehollandiae) co-occur with courtship dance wing flaps, while the ‘C’ type occurs when wings are still (Dalziell et al., 2013). In rats (Rattus norvegicus), 50 kHz ultrasonic vocalizations are produced during locomotor activity, while 20 kHz vocalizations occur when rats are immobile (Laplagne and Elías Costa, 2016). Similarly, bats in flight coordinate the production of echolocation sounds with particular phases of their wing beats (Suthers et al., 1972; Wong and Waters, 2001; Holderied and von Helversen, 2003). Thus, the adaptive coordination achieved during vocal development must also include other developing motor systems as important factors, notably those related to posture and locomotion. How this ‘vocal-locomotor’ coordination is accomplished over the course of development is not well understood.

In humans, current theoretical frameworks posit that the different motor systems are interactive over the course of development; one system can influence another in ways that change as a function of time (Thelen, 1991; Adolph, 2008; Iverson, 2010; Adolph and Robinson, 2015; Libertus and Hauf, 2017). In support of this, there is evidence that locomotor skills at one time point predict speech ability months or even years later (LeBarton and Iverson, 2013; Walle and Campos, 2014; Wang et al., 2014; He et al., 2015; LeBarton and Iverson, 2016; Libertus and Violi, 2016; Walle, 2016; Garrido et al., 2017; Libertus and Hauf, 2017; Salavati et al., 2017; West et al., 2019). However, only a handful of empirical studies investigated human infant vocal-locomotor coordination (Ejiri and Masataka, 2001; Fagan and Iverson, 2007; Abney et al., 2014; Berger et al., 2017). These studies showed that infant production of pre-linguistic vocalization is highly sensitive to body movement. For example, infants who are beginning to transverse their environment (i.e., crawling) are less likely to vocalize during locomotion than while they are sitting (Ejiri and Masataka, 2001; Fagan and Iverson, 2007; Abney et al., 2014; Berger et al., 2017). What is missing is an understanding of how vocal and locomotor outputs are coordinated in real-time in infants and how this coordination may change over the course of development.

Key to understanding these developmental dynamics is to also identify physiological conditions that may promote (or potentially inhibit) coordination between different motor outputs. One candidate is the state of arousal, a product of the autonomic nervous system and relevant for a range of behaviors (Pfaff, 2006). An animal would be said to exhibit a high arousal state if it is more alert to sensory stimuli, more motorically active and more reactive (Pfaff, 2006). The role of arousal is essentially to allocate metabolic energy (i.e., to prepare the body for action). As it relates to vocal production, arousal modulates respiration, which in turn provides the power for vocal output. Humans, for example, exhibit an increase in arousal--as measured by heart rate--prior to speaking (Lynch et al., 1980; Linden, 1987). In developing individuals, variable and spontaneous behaviors are ubiquitous, providing the scaffolding for more complex and organized behaviors later in life (Blumberg et al., 2013). These early behaviors, including vocal output and other bodily movements, primarily reflect the interplay between the infants’ arousal states, sensorimotor coordination and biomechanical conditions (Robinson et al., 2000). Thus, investigating the relationship between arousal fluctuations and the development of vocal and locomotor behaviors may prove to be illuminating.

Using marmoset monkeys as a model, here we investigate the relationship between vocal and locomotor systems and arousal levels during postnatal development. In the vocal domain, infant marmosets spontaneously produce sequences of immature and mature vocalizations, and these are linked to real-time changes in arousal levels (Zhang and Ghazanfar, 2016). Over the course of approximately two months, infants exhibit changes in the acoustic properties of their vocalizations that reflect a transition from producing mostly immature-sounding contact calls (e.g., cries) to mature-sounding contact calls (e.g., phees) (Takahashi et al., 2015; Zhang and Ghazanfar, 2016; Teramoto et al., 2017; Zhang and Ghazanfar, 2018). As in humans (Goldstein and Ja, 2008), this transition is facilitated by, and dependent upon, social reinforcement from parents (Takahashi et al., 2015; Gultekin and Hage, 2017; Takahashi et al., 2017; Gultekin and Hage, 2018). Moreover, these parallels with human prelinguistic development occur in the same life history stage (early infancy) (Ghazanfar and Liao, 2018). In the postural and locomotor domains, marmoset monkey development transitions from immature to mature forms in a pattern that is also similar to human development (e.g., righting reflex before sitting, and crawling before walking) (Wang et al., 2014; Braun et al., 2015; Schultz-Darken et al., 2016).

We address three fundamental questions: (1) Does one motor system – vocal or postural-locomotor – mature first or do they follow an overlapping trajectory? (2) How are vocal-postural-locomotor systems coordinated and do these coordination dynamics shift across development? (3) How do real-time fluctuations in arousal relate to vocal production, locomotion and their coordination?

Results

During their first two months of postnatal life, we measured infant marmoset behavior in a controlled context for 10 minutes approximately every 2 days. In each session, individuals were placed in a testing box in an experiment room that was outside visual and auditory range of their family groups. This brief ‘isolation’ context is a standard testing paradigm used to elicit vocalizations (Takahashi et al., 2015; Zhang and Ghazanfar, 2016) and to study the postures and locomotion of marmoset infants (Wang et al., 2014; Braun et al., 2015). The subjects were seven marmosets (three females) from three different parental pairs (two sets of twins, one set of triplets). We recorded the behaviors of each subject across ~30 sessions (28–33 per subject for a total of 220 sessions).

For vocal behaviors, we focused on two types of contact calls – cries and phees (Figure 1A). As described previously (Takahashi et al., 2015), cries are immature contact calls that have a short duration and noisy spectral properties (i.e., high Wiener entropy); phees are mature-sounding contact calls that have a longer duration and tonal spectral properties (i.e., low Wiener entropy). Cries transform into phees over the course of development (Takahashi et al., 2015; Zhang and Ghazanfar, 2016; Takahashi et al., 2017).

Vocal, postural and locomotor behaviors in developing marmosets.

(A) Spectrograms of immature (cries) and mature (phees) contact calls produced by infant marmosets in an isolated social context. (B) Cartoons representing the five types of posture behaviors. (C) Cartoons representing the five types of locomotor behaviors.

https://doi.org/10.7554/eLife.41853.002

We measured five types of postural behaviors – righting reflex, head raising, forelimb support, hindlimb support, and hanging (Wang et al., 2014; Braun et al., 2015) (Figure 1B). The righting reflex is when infants re-establish their body orientation so that their hands and feet are on the ground; head raising is when infants lift their head off the ground and look forward or up; forelimb support is when infants sit on the ground with their hands touching the ground; hindlimb support is when infants sit on the ground with their hands off the ground; hanging is when infants grasp elements in their environment (e.g., bars of the testing box) so that their hands and feet do not touch the ground.

We measured five types of locomotor behaviors – crawling, digging, jumping, climbing, and walking (Wang et al., 2014; Braun et al., 2015) (Figure 1C). Crawling is when infants move forward on the ground with their stomach touching the ground; digging is when infants move their hands back and forth across the ground; jumping is when infants push themselves off the ground or cage to move from one location to another; climbing is when infants traverse across the cage; walking is when infants traverse across the ground in a standing orientation.

Finally, for all seven infants, we concurrently measured arousal levels by acquiring heart rates during the sessions using non-invasive surface electrocardiography (Borjon et al., 2016; Zhang and Ghazanfar, 2016).

The vocal system matures before postural and locomotor systems

By measuring both vocal and postural-locomotor behaviors longitudinally in developing marmosets, we determined how these motor systems changed relative to one another. We first classified each behavior as immature or mature by measuring how their use shifted across development. In the vocal domain, the proportion of time producing cries decreased across development, while the proportion of time producing phees increased (Figure 2A; Table 1; Appendix 1.1) (Takahashi et al., 2015; Zhang and Ghazanfar, 2016). As indicated by physiology and biomechanics (Takahashi et al., 2015; Zhang and Ghazanfar, 2016; Zhang and Ghazanfar, 2018), cries were categorized as an immature contact call, and phees were categorized as a mature version of the contact call. The proportion of postural time engaged in righting reflexes decreased across development, while the proportion of postural time engaged in hindlimb support increased (Figure 2A; Table 1; Appendix 1.2). Thus, righting reflexes were categorized as an immature posture, and hindlimb support was categorized as a mature posture. The proportion of locomotor time engaged in crawling decreased across development, while the proportion of locomotor time engaged in walking increased (Figure 2A; Table 1; Appendix 1.3). Thus, crawling was categorized as an immature locomotor behavior, while walking was categorized as mature. These postural-locomotor classifications for developing marmosets are consistent with previous findings (Wang et al., 2014).

Contact calls mature before postural and locomotor behaviors.

(A) Changes in the proportion of time spent in specific vocal, postural, and locomotor behaviors across development. Behaviors that increase or decrease are denoted with asterisks to represent significance values. (B) Developmental trajectories of contact call maturation (i.e., mature phee calls relative to immature cries), postural maturation (i.e., hindlimb support relative to righting reflex), and locomotor maturation (i.e., walking relative to crawling). Points represent session values, gray curves represent the cubic spline fits for individual marmosets, and black curves represent the population cubic spline fits. (C) Illustration of hypothesized sequences between contact call and postural/locomotor developmental trajectories and (D) an overview of the observed developmental trajectories. (E) Comparison of the different developmental time courses. Bars and whiskers represent mean ±2 SE. Significance values are represented by ‘*' (p<0.05), ‘**' (p<0.01), ‘***' (p<0.001) and ‘****' (p<0.0001).

https://doi.org/10.7554/eLife.41853.003
Table 1
Results of linear mixed models (LMMs) to test whether proportion of time spent in vocal-postural-locomotor behaviors changes with postnatal day.

For each model, the proportion of vocal, postural, or locomotor time (per postnatal day) spent engaged in a behavior is the dependent variable, postnatal day is the fixed effect, and infant identity is the random effect. For each behavior category, a Bonferroni-Holmes correction was applied to adjust p-values.

https://doi.org/10.7554/eLife.41853.004
Behaviorβ (SE)T valueAdjusted PClassification
Vocal behaviors (n = 192 observation days)
Cry−0.0128 (0.0019)6.880.0011immature
Phee0.0129 (0.0017)7.430.0004mature
Postural behaviors (n = 201 observation days)
Forelimb support−0.0010 (0.0020)0.51>0.05NA
Hanging−0.0018 (0.0011)1.67>0.05NA
Hindlimb support0.0130 (0.0034)3.780.0373mature
Raising head−0.0021 (0.0009)2.32>0.05NA
Righting reflex−0.0081 (0.0015)5.540.0020immature
Locomotor behaviors (n = 191 observation days)
Climbing−0.0023 (0.0011)2.10>0.05NA
Crawling−0.0160 (0.0019)8.350.0019immature
Digging0.0033 (0.0011)3.03>0.05NA
Jumping0.0005 (0.0002)2.37>0.05NA
Walking0.0147 (0.0012)11.95<0.0001mature

We used the immature-mature classifications to summarize developmental changes using a ‘maturity index’ (see Materials and methods)--a value that represents the proportion of mature behaviors observed relative to all immature and mature behaviors across postnatal days. Values below 0.5 indicate that immature behavior is more common, and values above 0.5 indicate that mature behavior is more common. The population-level developmental trajectories of all three motor systems (vocal = 192 observation days, postural = 201 observation days, locomotor = 191 observation days) followed s-shape patterns that ended in maturity indices around one, which represents adult-like motor outputs (Figure 2B).

We tested whether the developmental time courses of the vocal system and postural-locomotor systems overlap (Figure 2C). The null hypothesis is that vocal and postural-locomotor systems transition from immature to mature forms around the same time. This would suggest that the development of the different motor systems reflect a global process of neural or physiological maturation (Gesell, 1929; McGraw, 1943). The alternative hypothesis is that the postural-locomotor system develops either before or after the vocal system. If postural-locomotor behaviors were to mature first, it would suggest that a developed post-cranial body is a prerequisite for producing mature contact calls.

Supporting the ‘vocal behavior develops first’ hypothesis, contact calls transitioned to a more mature form around 10 postnatal days, while postural and locomotor behavior transitioned around postnatal days 19 and 21, respectively (Figure 2D). A generalized linear mixed model (GLMM) showed that the relationship between maturity indices and postnatal day fits a logistic regression, with maturity indices increasing with age (n = 519 observation days, β ±SE = 0.21±0.03, z = 6.33, p<0.0001; Appendix 1.4). The same model showed that the vocal maturity indices were larger than the maturity indices of postural (β ±SE = −2.03±0.58, z = 3.48, p=0.0005; Appendix 1.4) and locomotor behavior (β ±SE = −2.86±0.98, z = 2.93, p=0.0034; Appendix 1.4) (Figure 2E). In other words, vocal development occurred earlier than postural and locomotor development. However, simply because vocal-postural-locomotor systems follow different trajectories does not mean that they do not interact in real-time. The question of whether vocal-locomotor coordination changes over the course of development is addressed next.

Mature contact call production and locomotor activity become increasingly coordinated during development

Infant marmosets require considerable muscular effort to produce mature contact calls (phees) (Takahashi et al., 2015; Zhang and Ghazanfar, 2016; Zhang and Ghazanfar, 2018), and adult marmosets tend to produce mature contact calls during periods of reduced locomotor activity (Borjon et al., 2016). Therefore, our overarching hypothesis was that locomotor activity inhibits mature contact call production. Using linear mixed effect models (LMMs), we found initial support for this hypothesis when examining the relationships between vocal acoustic parameters, locomotor activity and postnatal day (Figure 3A,B). Call duration was negatively associated with locomotor activity (n = 9609 calls, β ±SE = −0.56±0.13, t = 4.25, Bonferroni-Holm adjusted p=0.0001; Appendix 1.5) while Wiener entropy (higher entropy means noisier) was positively associated with locomotor activity (n = 9606 calls, β ±SE = 0.93±0.12, t = 7.68, Bonferroni-Holm adjusted p<0.0001; Appendix 1.6). In other words, on average, infant marmosets produced short, noisier calls when locomotor activity was high, and longer, more tonal adult-like calls when locomotor activity was reduced.

Figure 3 with 1 supplement see all
Mature contact calls and locomotor activity become increasingly coordinated across development.

(A) Plots of contact call acoustic parameters and corresponding locomotor activity in one infant. Green pluses represent immature contact calls (cries), and purple circles represent mature contact calls (phees). (B) Plots of contact call acoustic parameters and locomotor activity levels during contact calls. Points represent subject averages (one color per subject) per day and fit with a linear regression (black dotted line) with 95% confidence intervals (red dotted lines). (C) Locomotor activity from 20 s before to 20 s after contact call production. (C, top) Subject averages (one line per subject) of cubic splines fit to locomotor activity fluctuations surrounding cries (green lines) and phees (purple lines). (C, bottom) Colored cubic splines fit to population data are plotted with a 95% bootstrapped confidence interval (green – cries, purple – phees). Red line indicates observed values outside of the 95% threshold of the bootstrapped significance test. (D) Fluctuations in locomotor activity during contact calls across development. (D, top) Locomotor activity during cries (green pluses) and phees (purple circles) for all observation sessions. (D, bottom) Colored cubic splines fit to population data are plotted with a 95% bootstrapped confidence interval (green – cries, purple – phees). A red line indicates observed values outside of the 95% threshold of the bootstrapped significance test, and a black line indicates observed values inside the 95% threshold.

https://doi.org/10.7554/eLife.41853.005

Next, we tested three hypotheses about the developmental dynamics of vocal-locomotor coordination. We know that by 1–2 months of age, marmosets only produce mature sounding contact calls (Takahashi et al., 2015; Zhang and Ghazanfar, 2016; Zhang and Ghazanfar, 2018). One hypothesis is that they simply learn to stop moving when they need to produce a contact call. Doing so would eschew any potential physiological constraints on the production of mature sounding contact calls. In such a scenario, we would predict the number of contact calls produced during movement would decrease, while those produced during periods of immobility would increase. An alternative hypothesis is that as marmosets grow bigger, they become more capable of producing mature sounding contact calls while moving. This outcome would suggest that infants overcome any potential physiological constraints on vocal-locomotor coordination. In such a scenario, we would predict that the number of contact calls produced during movement would increase, while those produced during periods of immobility would decrease.

To test the above hypotheses, we first needed to estimate locomotor activity on a continuous scale. These estimates were extracted from frame-by-frame pixel differences in the video footage of the sessions (Figure 3—figure supplement 1). We first summarized the temporal real-time dynamics of locomotor activity surrounding vocal production events (Figure 3C). We found that locomotor activity started to increase above the 95% threshold of the bootstrap significance test ~5 s before the immature call onsets and then decrease back to inside the 95% threshold ~10 s after call offsets. In contrast, locomotor activity started to decrease below the 95% threshold of the bootstrap significance test ~8 s before mature contact call onsets and then increase back to inside the 95% threshold ~10 s after call offsets. A LMM confirmed that the average locomotor activity was higher during immature contact calls as compared to mature contact calls (n = 9609 calls, β ±SE = −0.07±0.01, t = 5.19, p<0.0001; Appendix 1.7). When these temporal dynamics were mapped across postnatal days, we found that locomotor activity during immature contact call production remained around or above the randomized expected levels of locomotor activity. In contrast, early in postnatal life, mature contact call production occurred when locomotor activity was below the 95% threshold of the bootstrap significance test; however, their production gradually became coordinated with increased levels of locomotor activity (Figure 3D). Locomotor activity during mature contact calls was higher during late development (PND 52–61) as compared to early development (PND 1–10) (Early = 914 calls; Late = 932 calls; β ±SE = 0.04±0.01, t = 3.05, p=0.0056; Appendix 1.8).

These data therefore support the alternative hypothesis, that the potential constraints of producing mature contact calls during movements are mitigated as the infants get older. Thus, even though postural-locomotor maturity does not appear to be a prerequisite for vocal development (Figure 2), vocal-locomotor coordination is still an important component of marmoset monkey motor development. This finding then raises the question of how real-time fluctuations in physiological condition may predict call production and locomotor activity in developing infants. Presumably, infants in elevated states of arousal are more likely to engage in these motor behaviors. The question of whether temporal associations between motor output and arousal levels change over the course of development is addressed next.

Mature call production and locomotion occurs during elevated arousal levels

We first tested hypotheses about the developmental dynamics of arousal state during marmoset contact call production. One very basic hypothesis is that the production of mature contact calls requires elevated arousal levels more than does the production of immature contact calls (Teramoto et al., 2017). Being in an elevated arousal state also means that individuals may have more respiratory power needed to generate mature sounding calls (Borjon et al., 2016; Zhang and Ghazanfar, 2016). Another hypothesis (that doesn’t preclude the first one) is that mature contact calls are more likely to occur during elevated arousal levels earlier in development. This would be consistent with computational models that suggest that younger infants are unable to effectively coordinate the elements of their vocal apparatus to produce mature contact calls and may require enhanced respiratory power to do so (Takahashi et al., 2015; Teramoto et al., 2017); indeed, there is now empirical support for the link between respiratory power and mature contact call production (Zhang and Ghazanfar, 2018). Conversely, a third hypothesis is that mature contact calls are more likely to occur during elevated arousal states later in development than earlier. This could occur if, as the infants get older, their overall motivation changes (i.e., less stressed by isolation) and so higher arousal levels are needed to motivate vocal production. An elevated arousal state may also enable respiratory power and laryngeal tension (Zhang and Ghazanfar, 2016); this could lead to better coordination between vocal production and locomotion (Figure 3C,D).

We found that mature contact calls occurred during elevated levels of arousal, a pattern that become more pronounced as development proceeded. A summary of the real-time dynamics of heart rate fluctuations in the 10 s before to 10 s after vocal production revealed that marmosets produced immature contact calls when heart rate percentiles dropped below the 95% threshold of the bootstrap significance test (Figure 4A). Production of mature contact calls, on the other hand, occurred when heart rate percentiles went above the 95% threshold of the bootstrap significance test. In both cases, changes in arousal levels were ‘global’, meaning that the change in arousal occurred well before (at least 10 s) the start of the call. A LMM confirmed that heart rate percentiles were higher during mature contact calls than during immature contact calls (n = 6215 calls, β ±SE = 3.38±1.45, t = 2.32, p=0.0276; Appendix 1.9). The developmental dynamics of arousal levels during vocal production further supported this main effect (Figure 4B). Early in postnatal life, heart rate percentiles during immature contact calls were inside the 95% threshold of the bootstrap significance test, but at around two weeks, heart rates during immature contact calls started to decrease below the 95% threshold. In contrast, mature contact calls were produced when heart rate percentiles were at, or above, the 95% threshold of the bootstrap significance test during the first two months of postnatal life. There was a marginal increase in heart rate percentiles during mature contact calls during late development (PND 52–61) as compared to early development (PND 1–10) (Early = 490 calls; Late = 780 calls; β ±SE = 6.28±3.48, t = 1.81, p=0.0795; Appendix 1.10).

Figure 4 with 1 supplement see all
Mature contact calls and locomotor activity are more likely to occur during elevated arousal levels later in development.

(A) Arousal fluctuations from 10 s before to 10 s after contact call production. (A, top) Subject averages (one line per subject) of cubic splines fit to arousal level fluctuations surrounding cries (green lines) and phees (purple lines). (A, bottom) Colored cubic splines fit to population data are plotted with a 95% bootstrapped confidence interval (green – cries, purple – phees). Red line indicates observed values outside of the 95% threshold of the bootstrapped significance test. (B) Fluctuations in arousal levels during contact calls across development. (B, top) Arousal levels during cries (green pluses) and phees (purple circles) for all observation sessions. (B, bottom) Colored cubic splines fit to population data are plotted with a 95% bootstrapped confidence interval (green – cries, purple – phees). Red line indicates observed values outside of the 95% threshold of the bootstrapped significance test. (C) Arousal fluctuations from 10 s before to 10 s after locomotor activity. (C, top) Subject averages (one line per subject) of cubic splines fit to arousal level fluctuations surrounding locomotor events (orange lines). (C, bottom) A colored cubic spline (orange) fit to population data is plotted with a 95% bootstrapped confidence interval. Red line indicates observed values outside of the 95% threshold of the bootstrapped significance test. (D) Fluctuations in arousal levels during locomotor activity across development. (D, top) Arousal levels during locomotor events (orange circles) for all observation sessions. (D, bottom) A colored cubic spline (orange) fit to population data is plotted with a 95% bootstrapped confidence interval. A red line indicates observed values outside of the 95% threshold of the bootstrapped significance test, and a black line indicates observed values inside of the 95% threshold.

https://doi.org/10.7554/eLife.41853.007

We next tested hypotheses about the developmental dynamics of arousal during locomotor activity. One basic hypothesis is that locomotion occurs during elevated arousal states, which is a pattern well supported by studies of heart rate during physically demanding tasks (Rotstein et al., 2004; Baker et al., 2008). Another hypothesis (that is not mutually exclusive of the first one) is that locomotor activity is more likely to occur during elevated arousal levels earlier in development. One interpretation of this result is that, as with vocal production, younger infants are less able to coordinate their body movements and require enhanced respiratory power to do so. Conversely, a third hypothesis is that locomotor activity is more likely to occur during elevated arousal levels later in development than earlier. As with vocal production, this could occur if, as the infants get older, their overall motivation changes (i.e., less stressed by isolation) and so higher arousal levels are needed to motivate movement. An elevated arousal state may also enable respiratory power and laryngeal tension to enable better coordination between vocal production and locomotion (Figure 3C,D).

Like vocal production, we found that locomotor activity occurred during elevated arousal levels, a pattern that appeared to become more pronounced later in development. A summary of the real-time dynamics of heart rate fluctuations in the 10 s before to 10 s after locomotion events revealed that these events occurred when heart rate percentiles were elevated above the 95% threshold of the bootstrap significance test (Figure 4C). As with vocal production, this change in arousal was ‘global’, meaning that the change in arousal happened well before (at least 10 s) the start of movement. The developmental dynamics of arousal levels during locomotor activity further supported this main effect result (Figure 4D). Early in postnatal life, heart rate percentiles during locomotor activity were within the 95% threshold of the bootstrap significance test, but continued to increase throughout the first two months of postnatal life. There was a marginal increase in heart rate percentiles during locomotor events during late development (PND 52–61) as compared to early development (PND 1–10) (Early = 752 calls; Late = 566 calls; β ±SE = 8.12±4.04, t = 2.01, p=0.0580; Appendix 1.11).

Our data suggest that mature contact call production and locomotor activity are both associated with elevated levels of arousal, an association that becomes more pronounced with age. As such, arousal state may be an important predictor of whether infant marmosets coordinate these different motor outputs. We tested this idea next.

Coordination of mature contact call production with locomotor activity occurs during elevated levels of arousal

Over the course of infant marmoset development, both locomotor activity and arousal levels predicted whether an infant produces a mature contact call over an immature call. Low levels of locomotor activity support mature contact call production early in development, and elevated arousal levels support mature contact call production later in development. The missing piece of the puzzle is whether there is a connection between vocal-locomotor coordination and arousal state. Here, we test, and find support for, the hypothesis that arousal levels during mature contact call production are elevated when marmosets are moving around (Figure 5A). Locomotor activity during mature contact call production was positively associated with heart rate percentiles (n = 3966 calls, β ±SE = 9.66±4.05, t = 2.39, p=0.0208; Appendix 1.12; Figure 5B). In other words, this result suggests that individuals in an elevated arousal state are better able to coordinate mature vocal production with locomotion. Specifically, this positive association characterized infants that were one month old (1–30 days; n = 1705 calls, β ±SE = 16.32±6.60, t = 2.47, Bonferroni-Holm adjusted p=0.0382; Appendix 1.12), but not two months old (31–61 days; n = 2261 calls, β ±SE = 2.63±5.24, t = 0.50, Bonferroni-Holm adjusted p=0.6220; Appendix 1.12). This means that a positive association between locomotor activity and heart rate percentiles was not simply a by-product of age. Instead, there may be a particularly robust relationship between locomotor activity and heart rate during early infanthood when individuals are transitioning from producing immature to mature contact calls.

Coordination of mature contact calls with locomotor activity occurs during elevated arousal levels.

(A) Plots of locomotor activity and heart rate percentiles during contact calls in one infant. Purple circles represent mature contact calls (phees). (B) Plot of heart rate percentiles and locomotor activity during mature contact calls (phees). Points represent subject averages (one color per subject) per day and fit with a linear regression (black dotted line) with 95% confidence intervals (red dotted lines).

https://doi.org/10.7554/eLife.41853.009

Discussion

Vocal development is a dynamic process that involves the interaction of multiple systems, from the biomechanics and muscles of the vocal apparatus to the nervous system and the social environment (Thelen, 1991; Teramoto et al., 2017). Using developing marmoset monkeys, we sought to understand how these processes related to vocal output are influenced by other systems of motor behavior, specifically posture and locomotion. First, we examined when the transition from immature to mature calls occurs relative to the transition from immature to mature body postures and locomotion. Second, we examined the putative temporal coordination of vocal production with locomotor activity, and whether such coordination changes over the course of development. Finally, we investigated whether fluctuating arousal levels (estimated from heart rate) predicts vocal and locomotor output. We found that marmoset monkey vocalizations develop sooner than postural-locomotor skills, that locomotor activity gradually becomes coordinated with the production of mature-sounding contact calls, and that this vocal-locomotor coordination occurs during elevated levels of arousal.

Head-to-tail sequence of development

The development of the vocal and postural-locomotor systems in marmoset monkeys, at a first approximation, seems to follow similar trajectories (Wang et al., 2014; Takahashi et al., 2015). By investigating both motor systems in the same individuals longitudinally, however, we showed that vocal behavior matures sooner than either posture or locomotion. Marmosets transitioned to producing a higher proportion of mature-sounding contact calls (phees) than immature-sounding contact calls (cries) around postnatal day 10, whereas posture and locomotion transitioned to more mature forms around days 18 and 21, respectively. This finding implies that the ability to produce adult-like sounds is not contingent upon advanced postural and locomotor control. These findings also parallel the developmental sequence of human prelinguistic vocal output, posture, and locomotion. Human infants begin producing babbling sounds (consonant-vowel combinations) at around 2 to 4 months (Vihman, 2014), unsupported upright sitting around 4 to 5 months (Adolph, 2008), and walking without support around 11 months (Adolph, 2008). Our data suggest that, like humans and other animals (Starck and Ricklefs, 1998; Muir, 2000; Adolph, 2008), marmoset monkey motor development takes the form of a head-to-tail sequence.

Why do marmoset monkey and human infants both transition to producing more adult-like sounds before they transition to sitting in upright postures and walking? One explanation is that an initial investment in the vocal system allows infants to elicit caregiver attention (e.g., carrying and food sharing) more effectively, thereby prolonging the amount of time they need to develop their locomotor autonomy. This explanation makes sense given the developmental strategy (i.e., altricial) and social system (i.e., cooperative breeding) of humans and marmosets (Snowdon, 1996; Hrdy, 2007; Solomon and French, 2007; Burkart et al., 2009; Lukas and Clutton-Brock, 2012). For altricial species like marmoset monkeys (relative to other nonhuman primates), early development is an energetically costly period because individuals are unable to fully regulate their body temperature, and yet, altricial infants must invest energy to grow at a high rate and refine their motor skills (Rosenblatt, 1976; Case, 1978; Derrickson, 1992; Blumberg and Sokoloff, 1998; Starck and Ricklefs, 1998; Muir, 2000; Blumberg, 2001; Schilling, 2005; Price and Dzialowski, 2018). In marmosets and humans, locomotor and physiological constraints (such as control of arousal levels) may be overcome by receiving care and physical contact from both maternal and non-maternal adults (Case, 1978; Snowdon, 1996; Hrdy, 2007); infants elicit such contact and care by producing vocalizations (Locke, 2006; Zuberbühler, 2012; Ghazanfar and Takahashi, 2017). And yet, not all vocalizations are created equal. Previous work in marmosets and humans suggests that more adult-like sounds elicit caregiver attention better than do immature sounds (Gros-Louis et al., 2006; Takahashi et al., 2016). By investing first in developing mature-sounding contact calls, infants may be ‘buying the time’ they need to learn how to move about independently in their environment.

Locomotion as a potential constraint on vocal development

The ability to coordinate biomechanical features of the vocal system – breathing, thoracic pressure and vocal fold tension – is critical for producing species-typical vocalizations (MacLarnon and Hewitt, 1999; Maclarnon and Hewitt, 2004; Takahashi et al., 2015; Zhang and Ghazanfar, 2016; Teramoto et al., 2017). In developing marmosets, a mismatch between biomechanical dynamics of the vocal system is thought to generate immature cries instead of mature contact calls (Teramoto et al., 2017). Our study indicates that locomotor activity is one force that can potentially disrupt vocal production. Higher levels of body movements co-occurred with the production of immature-sounding contact calls (i.e., those with a short duration and high entropy), while lower levels of movement co-occurred with mature-sounding contact calls with long durations and lower entropy. These results are consistent with what is known about how movement affects respiration. Vigorous motor activities, like running, speed up breathing cycles (Wasserman et al., 1973; Bramble and Carrier, 1983), resulting in articulation deficits (Sundberg et al., 1991; Price et al., 2006; Baker et al., 2008; Orlikoff, 2008). The finding that calls produced during movement were shorter and noisier indicates that very young marmosets lacked adequate respiratory power (Zhang and Ghazanfar, 2018).

Arousal levels during vocal-locomotor coordination

A unique aspect of our study design is that we could test how real-time fluctuations in arousal related to vocal production, locomotion, and the coordination of these different motor outputs. Infant marmosets produced mature-sounding contact calls during elevated arousal levels and immature ones during low arousal levels, a finding that is consistent with previous work on infant and adult marmosets (Borjon et al., 2016; Zhang and Ghazanfar, 2016). Similar to mature-sounding contact calls, infant marmosets also tended to engage in locomotion during elevated arousal states. Moreover, real-time variability in arousal state indicated whether mature-sounding contact calls co-occurred with locomotor activity. Infant marmosets exhibited elevated arousal levels when mature-sounding contact calls were produced during movement, and decreased arousal levels when such calls were produced during periods of immobility. These results are consistent with the hypothesis that elevated arousal may help to overcome physiological demands (e.g., respiratory power) of producing vocalizations while moving. Similar associations have been observed in humans engaged in physically demanding tasks (Rotstein et al., 2004; Baker et al., 2008). For example, heart rate increases at a faster rate during sustained exercise in adults who are engaged in a speech task than in adults engaged in a non-speech task (Baker et al., 2008). Our study design precludes testing causality but does suggest that arousal state is a key player in the coordination between vocal and locomotor systems during development.

The inhibition of motor activity during mature contact call production shifted over the course of development. That is, by the time that infant marmoset monkeys stopped producing immature-sounding calls, they no longer showed decreased movement during mature contact call productions. This developmental shift suggests that marmosets gradually improved their ability to coordinate locomotor behaviors with vocal production. From a biomechanical perspective, this improvement suggests that marmosets acquired the ability to better control the vocal apparatus during movement. From a neural perspective, this improvement suggests that more experience with vocal behaviors and/or locomotion leads to better coordination between these motor systems. In either case, arousal state appears to be the common currency by which vocal-motor coordination emerges. One intriguing possibility is that these shifting dynamics of the autonomic nervous system create the scaffolding by which mature social behavior can emerge (Porges and Furman, 2011).

Despite being identified as a critical line of inquiry over 30 years ago (Tipps et al., 1981; Yingling, 1981), to date, only a handful of empirical studies investigated infant vocal-locomotor coordination (Ejiri and Masataka, 2001; Fagan and Iverson, 2007; Abney et al., 2014; Berger et al., 2017). Our study using marmoset monkeys as a model for developmental processes represents one of the first to integrate longitudinal and second-by-second timescales to investigate vocal development from a ‘whole-body’ perspective. We believe that this timescale integration is key to characterizing the fine-grained dynamics that dictate how mature vocalizations emerge. Trade-offs in vocal-locomotor coordination is a potential dynamic that needs to be reckoned with as individuals grow up. By concurrently measuring heart rate, we show that processes related to autonomic arousal may enable individuals to cope with this trade-off.

Materials and methods

Subjects and housing

Request a detailed protocol

The subjects used in this study were seven (three females) infant common marmosets (Callithrix jacchus) from three different parental pairs (two sets of twins and one set of triplets,<2 months old). Subjects were born in captivity and lived with their family groups (mother, father and siblings). The colony room was maintained at a temperature of approximately 27° C with 50–60% relative humidity, and a 12:12 hr light-dark cycle. All subjects had ad libitum access to water and were supplied daily with standard commercial chow supplemented with fruit and vegetables.

Experimental design

Request a detailed protocol

The experimental protocol follow methods previously described (Borjon et al., 2016; Zhang and Ghazanfar, 2016). Infant marmosets were separated from their parents and placed in a testing box in an experiment room. The triangular, prism-shaped testing boxes were made of Plexiglas and wire (0.30 m x 0.30 m x 0.35 m). All observation sessions were conducted during daylight hours between postnatal day 1 to 61, and each observation session lasted 10 minutes. Subjects participated in a total of 220 observation sessions across the first two months of life (Subjects 1–7: 29, 29, 34, 34, 31, 31, 32). All experimental procedures were performed in compliance with the guidelines of the Princeton University Institutional Animal Care and Use Committee.

Vocal behavior data collection

Request a detailed protocol

Undirected vocalizations (i.e., socially isolated) were recorded using a Sennheiser MKH416-P48 microphone suspended 0.9 m above the testing box. The microphone signal was sent to a Mackie 402-VLZ3 line mixer whose output was relayed to a Plexon Omniplex and PC computer. We used the same custom MATLAB software established in previous research for computationally defining and segmenting infant marmoset vocalizations (Takahashi et al., 2015; Zhang and Ghazanfar, 2016). A researcher manually verified which calls were one of two types of contact calls – cries and phees (Figure 1A). As described previously (Takahashi et al., 2015), cries are contact calls that have a short duration and noisy spectral properties (i.e., high Wiener entropy); phees are contact calls that have a longer duration and tonal spectral properties (i.e., low Wiener entropy). We recorded audio for 192 of the 220 observation sessions (Subjects 1–7: 28, 29, 21, 20, 31, 31, 32) for a total of 10,956 calls (Subjects 1–7: 1,100, 801, 1,604, 1,021, 1,657, 2,199, 2,574). A custom-made MATLAB routine calculated two main acoustic properties for each call, the duration and Wiener entropy (Takahashi et al., 2016). Wiener entropy is a non-positive number that is calculated by taking the logarithm of the ratio between the geometric and arithmetic means of the values of the power spectrum for different frequencies (Tchernichovski et al., 2000). Wiener entropy represents the broadband properties of a signal’s power spectrum in which the closer the signal is to white noise, the higher (closer to zero) the entropy value.

Postural and locomotor behavior data collection

Request a detailed protocol

Postural-locomotor behavior was video recorded at 30 fps with a Plexon Cineplex. We recorded video for 215 of the 220 observation sessions and identified a total of 3195 instances of specific postural-locomotor behaviors (Subjects 1–7: 662, 635, 414, 574, 207, 427, 276). We manually scored the recorded videos to identify the onset and offset of behaviors using BORIS, an open-source event-logging software for video coding (Friard and Gamba, 2016). As the video frame rate is at 30 Hz, the onset and offset of behaviors had a maximum resolution of 1/30 s. Our definitions of postural-locomotor behaviors were based on prior literature (Wang et al., 2014). During video coding, frames without a clear posture or locomotion described in our ethogram were not assigned a behavior.

Postural behaviors were defined as instances that infants spent repositioning itself: forelimb support, hanging, hindlimb support, raising head, and righting reflex (Figure 1B). The righting reflex is when infants re-establish their body orientation so that their hands and feet are on the ground; head raising is when infants lift their head off the ground and look forward or up; forelimb support is when infants sit on the ground with their hands touching the ground (or cage); hindlimb support is when infants sit on the ground with their hands off the ground (or cage); hanging is when infants grasp the cage so that their hands and feet do not touch the ground. Locomotor behaviors were defined as instances when infants traversed the testing box: crawling, digging, jumping, climbing, and walking (Figure 1C). Crawling is when infants move forward on the ground with their stomach touching the ground; digging is when infants move their hands back and forth across the ground; jumping is when infants push themselves off the ground or cage to move from one location to another; climbing is when infants traverse across the cage; walking is when infants traverse across the ground in a standing orientation.

Continuous variability in locomotor activity (i.e., body movement) was assessed by investigating pixel differences in the video data (Figure 3—figure supplement 1), following methods used in earlier research (Borjon et al., 2016). Each video recording was split into segments of 30 frames (each 640 vs 400 pixels). We took the absolute difference of RGB values between the first and last frame of every second and divided by the total number of pixels. This value corresponded to the average difference in luminescence per pixel per second. A higher value indicates more pixel difference, signifying movement. Because absolute levels of locomotion differ across individuals and ages (e.g., movement of larger individuals causes larger pixel differences), it was necessary to re-scale locomotor activity levels. To do this, we converted all 1 s movement values to binary values with a 90th percentile threshold. Then, we use csaps in MATLAB to fit a cubic spline (smoothing parameter of 0.10) to the binary values in each 10 min observation session. In other words, locomotor activity ranged between zero (immobile) and one (mobile) so that comparisons could be made across individuals and ages.

Electrocardiography data collection

Request a detailed protocol

To quantify arousal fluctuations, we recorded heart rate for 149 of the 220 observation sessions (Subjects 1–7: 16, 14, 21, 19, 26, 25, 28). To record electrocardiographic (ECG) signal, we used two pairs of Ag-AgCl surface electrodes (Grass Technology) (Figure 4—figure supplement 1). Tethered electrodes were sewn into a soft elastic band, which was clasped around the animal’s thorax. One pair of electrodes was positioned on the dorsal thorax, and the other pair was positioned on the ventral thorax. We applied ECL gel on the surface of each electrode to improved signal-to-noise ratio. Infants were shaved around the thorax if needed. Each pair of electrodes was differentially amplified (x250) with the resulting signal sent to the Plexon Omniplex, where it was digitized at 40 kHz and sent to a personal computer for data acquisition. The strength of heart rate signals varies throughout observation sessions as animal movement alters the positioning of surface electrodes. As such, we chose the channel with the largest signal-to-noise ratio on a session-by-session basis. We manually identified and isolated motion artifacts or signal cutoffs. To minimize bias, we did the following for each observation session: signals from both ECG channels were divided into 10 s segments and signal pairs were presented in random order for visual inspection. Regions exhibiting signal loss were replaced with NaNs.

Following the visual screening, we down-sampled data from 40 kHz to 1,500 Hz to extract the cardiac signals. These signals were high-pass filtered at 15 Hz to preserve the rapid waveform of the heartbeat. The resulting signal was notch filtered at 60 Hz. Heart beats were detected using an adaptive threshold of 1 s duration to find cardiac spikes greater than the 95th percentile of the amplitude at each second of the signal. Occasionally, an artificial spike close to the actual heartbeat would be detected or a heartbeat would be missed. As such, we set inter-spike interval thresholds of 100 ms (600 beats/min) to 400 ms (150 beats/min), which are thresholds used in a previous study of marmoset autonomic activity (Borjon et al., 2016). If an inter-spike interval was less than 100 ms, the two spikes were substituted with a single spike located at the midpoint. If an inter-spike interval exceeded 400 ms, we replaced the signal with a NaN. To calculate heart rate, we constructed a binary series of heartbeat counts and convolved the resulting series with a 1 s Gaussian window. We only used heart rate data from observation sessions during which heart rate could be detected at least 50% of the time. Because heart rate can differ across individuals and ages, it was necessary to re-scale heart rate fluctuations so that comparisons could be made across individuals and ages. To do so, we converted all 1 s heart rate values to percentiles for each 10 min observation session (i.e., heart beat fluctuations were centered around the 50th percentile). In other words, heart rate percentile values ranged from zero (lowest heart rate level) to 100 (highest heart rate level) within each 10 min observation session.

Data analysis

Request a detailed protocol

All analyses were carried out in MATLAB (version R2019a) and R (version 3.5.3). To determine which motor behaviors were categorized as ‘immature’ and ‘mature’, we calculated the proportion of vocal, postural, or locomotor time engaged in specific behaviors. We used a series of linear mixed effect models (LMM, ‘lmer’ of the R package ‘lme4’; Bates et al., 2015) to test whether the proportion of time engaged in these behaviors changed based on postnatal day. We used the ‘lmerTest’ R package to determine the significance of the coefficients (Kuznetsova et al., 2017). In the LMMs to examine how vocal behavior changes across development, the dependent variable was the proportion of total vocal time (per daily observation session) engaged in a specific vocalization (cry or phee), the fixed effect was postnatal day, and the random effect was the infant subject (LMM equation: proportion of time ~day + (day|subject)). In the LMMs to examine how postural behavior changes across development, the dependent variable was the proportion of total postural time (per daily observation session) engaged in a specific posture (forelimb support, hanging, hindlimb support, raising head, or righting reflex), the fixed effect was postnatal day, and the random effect was infant subject (LMM equation: proportion of time ~day + (day|subject)). In the LMMs to examine how locomotor behavior changes across development, the dependent variable was the proportion of total locomotor time (per daily observation session) engaged in a specific posture (climbing, crawling, digging, jumping, or walking), the fixed effect was postnatal day, and the random effect was infant subject (LMM equation: proportion of time ~day + (day|subject)). We applied the Bonferroni-Holm method to correct for issues of multiplicity within each behavior type (vocal, postural, or locomotor behaviors), resulting in adjusted p-values with an alpha threshold level of 0.05. We report detailed outcomes of these regression models (e.g., model formulas, random effect variance, regression coefficients, standard errors, t-values, and p-values) in Appendix 1.1-1.3. Immature behaviors were those whose use decreased across development, and mature behaviors were those whose use increased across development. We used these immature-mature categories to calculate ‘maturity indices’ per session. This index ranged from 0 to 1 and was calculated as follows,

Maturity Index= mm+im

where m represents the percent of time spent engaged in mature behavior and im represents the percent of time engaged in immature behavior. A maturity index value less than 0.5 means that an individual produced more immature behavior, and a value greater than 0.5 means that an individual produced more mature behavior. Cubic splines (MATLAB csaps function) were fit to individual data (smoothing parameter of 0.03) and population data (smoothing parameter of 0.01) to determine the developmental trajectories of vocal, postural and locomotor maturity indices. Then, we used a logistic generalized linear mixed effect model (GLMM, ‘glmer’ of the R package ‘lme4’; Bates et al., 2015) to test whether maturation time courses differed between vocal and postural-locomotor behaviors. We used the ‘lmerTest’ R package to determine the significance of the coefficients (Kuznetsova et al., 2017). In this GLMM, the dependent variable was maturity index (per daily observation session, per behavior type—vocal, postural, locomotor), the fixed effect was postnatal day, and the random effect was infant subject (GLMM equation: maturity index ~behavior type+day + (behavior type|subject) + (day|subject)). In Appendix 1.4, we report detailed outcomes of this regression model (e.g., model formulas, random effect variance, regression coefficients, standard errors, t-values, and p-values), as well as the outcomes of this model per infant subject.

We sought to understand the real-time and developmental dynamics between vocal production and locomotor activity, as well as between arousal fluctuations and vocal-locomotor behavior. Then, to visualize the real-time dynamics of locomotion during call production, we extracted locomotor activity from −20 to 5 s surrounding call onsets and −5 to 20 s surrounding call offsets. To visualize the real-time dynamics of arousal fluctuations surrounding vocal-locomotor events, we extracted heart rate percentiles from −10 to 5 s surrounding call (or locomotor activity) onsets and −5 to 10 s surrounding call (or locomotor activity) offsets. To summarize the real-time dynamics for individual sessions, we fit a cubic spline to the session data (MATLAB csaps, smoothing parameter of 0.1), and then we fit a population spline to all session splines (smoothing parameter of 0.3). To visualize developmental dynamics, we used the average locomotor activity (or heart rate percentile) during each call (or locomotor activity event) to calculate the mean locomotor activity (or heart rate percentile) per individual session. Then, we fit cubic splines (smoothing parameter of 0.0001) to individual sessions to model the change in the population data across postnatal days.

To test the statistical significance of real-time and developmental dynamics, confidence intervals for the population splines were generated by randomly resampling (10 samples per session) from the signals used to generate individual session splines (real-time dynamics) or session averages (developmental dynamics). We fit a cubic spline to the resampled data to generate a population spline, and repeated the process 1000 times. The 95% confidence interval corresponds to the 2.5th and the 97.5th percentiles of the resampled population splines. To determine whether the real-time or developmental dynamics of locomotor activity and arousal fluctuations were significantly different from null expectations, we performed bootstrapped significance tests. For each session, we scrambled the order of call durations and inter-call interval durations. This allowed us to choose random segments equivalent in the number and length to the calls produced in that session while maintaining naturalistic spacing between the calls. We fit a cubic spline to session splines (for temporal analysis) or session averages (for developmental analysis) to generate bootstrapped population splines. We repeated this process 1000 times. The 95% threshold for significance corresponds to the 2.5th and 97.5th percentiles of the bootstrapped population splines. With this bootstrap procedure, we are taking into account data variability due to day, subject and call timing on the statistical estimates.

To complement our bootstrap procedure, we used a series of linear mixed effect models (LMM, ‘lmer’ of the R package ‘lme4’; Bates et al., 2015) to investigate associations between call production, locomotor activity, and arousal fluctuations. We used the ‘lmerTest’ R package to determine the significance of the coefficients (Kuznetsova et al., 2017). We used LMMs to investigate how average locomotor activity during a call predicted call duration and Wiener entropy. In these LMMs, the dependent variable was contact call acoustic parameter (call duration or Wiener entropy), the fixed effects was average locomotor activity (per call), and the random effect was infant subject nested within postnatal day (LMM equation: acoustic parameter ~locomotor activity + (locomotor activity|day/subject)). We applied the Bonferroni-Holmes method to correct for multiplicity issues associated with testing two acoustic parameters, resulting in adjusted p-values with an alpha threshold level of 0.05. In Appendix 1.5-1.6, we report detailed outcomes of these regression models (e.g., model formulas, random effect variance, regression coefficients, standard errors, t-values, and p-values), as well as the outcomes of these models per postnatal day and infant subject. To complement the real-time dynamic analyses, we used LMMs to compare locomotor activity and ANS activity between contact call types. For the LMM examining locomotor activity, the dependent variable was the average locomotor activity (per contact call), the fixed effect was contact call type (immature vs. mature), and the random effect was infant subject nested within postnatal day (LMM equation: locomotor activity ~call type + (call type|day/subject)). For the LMM examining ANS activity, the dependent variable was the average heart rate percentile (per contact call), the fixed effect was contact call type (immature vs. mature), and the random effect was infant subject nested within postnatal day (LMM equation: heart rate ~call type + (call type|day/subject)). In Appendices 1.7 and 1.9. we report detailed outcomes of these regression models for example model formulas, random effect variance, regression coefficients, standard errors, t-values, and p-values), as well as outcomes of the models per postnatal day and infant subject.

We also used LMMs to complement the developmental dynamics analyses. We ran an LMM to examine developmental changes in locomotor activity during mature contact calls, in which the dependent variable was the average locomotor activity (per contact call), the fixed effect was age group, and the random effect was infant subject nested within postnatal day (LMM equation: locomotor activity ~age group + (age group|day/subject)). Age group was split into early development (PND 1–10; infant contact calls transition to sounding more adult-like around PND 10) and late development (PND 52–61). In Appendix 1.8, we report detailed outcomes of this regression model, as well as outcomes for each infant subject separately. We ran an LMM to examine developmental changes in ANS activity during mature contact calls, in which the dependent variable was the average heart rate percentile (per contact call), the fixed effect was age group, and the random effect was infant subject nested within postnatal day (LMM equation: heart rate ~age group + (age group|day/subject)). In Appendix 1.10, we report detailed outcomes of this regression model, as well as outcomes for each infant subject separately. We also ran an LMM to examine developmental changes in ANS activity during locomotor activity events, in which the dependent variable was the average heart rate percentile (per locomotor event) the fixed effect was age group, and the random effect was infant subject nested within postnatal day (LMM equation: heart rate ~age group + (age group|day/subject)). In Appendix 1.11, we report detailed outcomes of this regression model, as well as outcomes for each infant subject separately.

Finally, we used a LMM to investigate how ANS activity during mature contact calls is associated with vocal-locomotor coordination. The dependent variable of this LMM was average heart rate percentile (per mature contact call), the fixed effect was average locomotor activity (per mature contact call), and the random effect was infant subject nested within postnatal day (LMM equation: heart rate ~locomotor activity + (locomotor activity|day/subject)). We also ran this model on two separate subsets of the data, the first month of development (postnatal day 1–30) and the second month of development (postnatal day 31–61). We applied a Bonferroni-Holmes method to correct for issues of multiplicity due to testing different time frames separately, resulting in adjusted p-values with an alpha threshold level of 0.05. In Appendix 1.12, we report detailed outcomes of these regression models (e.g., model formula, random effect variance, regression coefficients, standard errors, t-values, and p-values), as well as the outcomes for each postnatal day and subject separately.

Appendix 1

Linear Regression Models

We report below the model formulas, random effect variance, estimated regression coefficients, standard errors, test statistics and p values of the models reported in the main text.

1.1 Linear mixed models to predict the proportion of vocal time spent producing specific contact calls (cries or phees) from postnatal day. Results are associated with Figure 2A

1.1.1 Results for cries

Formula: Cry ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)8.116e-020.28489
Day8.116e-020.00446−1.00
Residual3.973e-020.19931

Number of obs: 192, groups: Subject, 7

Fixed effects:

EstimateStd. error dfT valuePr(>|t|)
(Intercept)0.5925730.1103384.6401885.3710.00377
Day−0.0128490.0018674.879883−6.8830.00109

1.1.2 Results for phees

Formula: Phee ~ Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)6.786e-020.26049
Day1.641e-050.00405−1.00
Residual3.992e-020.19981

Number of obs: 192, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.4071450.1013716.0407264.0160.006891
Day0.0128570.0017296.3378507.4340.000235

1.2 Linear mixed models to predict proportion of posture time spent engaged in specific postures (forelimb support, hanging, hindlimb support, raising head, or righting reflex) from postnatal day. Results are associated with Figure 2A

1.2.1 Results for forelimb support

Formula: ForelimbSupport ~ Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)2.555e-020.159854
Day1.414e-050.0037610.5
Residual1.250e-010.353497

Number of obs: 201, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.3773780.0726365.297525.1950.00294
Day−0.0010240.0020175.53592−0.5080.63123

1.2.2 Results for hanging

Formula: Hanging ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)1.008e-030.031752
Day1.030e-070.000321−1.00
Residual7.370e-020.271470

Number of obs: 201, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.1577810.03317210.8209464.7560.000621
Day−0.0018330.00109894.343478−1.6690.098491

1.2.3 Results for hindlimb support

Formula: HindlimbSupport ~ Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)9.179e-050.009581
Day7.557e-050.008693−1.00
Residual5.857e-020.242022

Number of obs: 201, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)−0.0366930.02784992.117003−1.3180.19091
Day0.012970.0034295.9420583.7820.00933

1.2.4 Results for raising head

Formula: RaiseHead ~ Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)3.447e-030.058707
Day1.221e-060.001105−1.00
Residual4.122e-020.203028

Number of obs: 201, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.1160940.0320547.0429293.6220.0084
Day−0.0021260.00091712.62578−2.3190.0379

1.2.5 Results for righting reflex

Formula: RightingReflex ~ Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)6.187e-030.07866
Day4.616e-060.002149−1.00
Residual9.023e-020.300387

Number of obs: 201, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.3865060.0453376.1605518.5250.000124
Day−0.0080660.0014568.71568−5.540.000405

1.3 Linear mixed models to predict proportion of locomotor time spent engaged in specific locomotor behaviors (climbing, crawling, digging, jumping, or walking) from postnatal day. Results are associated with Figure 2A

1.3.1 Results for climbing

Formula: Climb ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)4.207e-030.064858
Day3.389e-080.00018411.00
Residual7.232e-020.2689296

Number of obs: 191, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.1973060.0396510.0261264.9760.000552
Day−0.0023250.001109139.768818−2.0960.037849

1.3.2 Results for crawling

Formula: Crawl ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)3.671e-020.191607
Day1.313e-050.003624−0.83
Residual1.064e-010.326136

Number of obs: 191, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.849630.0817074.29213210.3980.000333
Day−0.0160480.0019214.816723−8.3550.000481

1.3.3 Results for digging

Formula: Dig ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)8.017e-040.028315
Day6.872e-060.002621−1.00
Residual1.135e-020.106525

Number of obs: 191, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)−0.0324130.0163478.552887−1.9830.0804
Day0.0032860.0010856.0502163.0280.0229

1.3.4 Results for jumping

Formula: Jump ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)1.708e-080.0001307
Day2.293e-070.0004789−0.30
Residual9.220e-040.0303643

Number of obs: 191, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)−0.00259600.00352118.830−0.7370.4619
Day0.00052260.00022069.1522.3690.0415

1.3.5 Results for walking

Formula: Walk ~Day + (Day | Subject)

Data: BehaviorProportions

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)3.903e-030.0624716
Day7.360e-070.00085790.98
Residual8.263e-020.2874487

Number of obs: 191, groups: Subject, 7

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)−0.0138530.04083311.37347−0.3390.741
Day0.0146830.00122918.15156411.9514.87e-10

1.4 Linear regression models to compare maturity indices between vocal and postural-locomotor behavior

1.4.1 Population-level logistic generalized linear mixed model results

Formula: Index ~ as.factor(BehaviorTypeDummy)+Day + (BehaviorTypeDummy | Subject) + (Day|Subject)

Data: MaturityIndices

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject(Intercept)0.7607310.87220
BehaviorTypeDummy1.3130621.14589−1.00
Subject.1(Intercept)0.0221070.14868
Day0.0039830.063110.99

Number of obs: 519, groups: Subject, 7

Fixed effects:

EstimateStd. errorZ valuePr(>|z|)
(Intercept)−1.921290.45544−4.2192.46e-05
as.factor(BehaviorTypeDummy)1−2.033260.58452−3.4790.000504
as.factor(BehaviorTypeDummy)2−2.863050.97803−2.9270.003419
Day0.208780.032996.3292.47e-10

1.4.2 Logistic generalized linear model results per Subject

Formula: Index ~ as.factor(BehaviorTypeDummy)+Day

Data: MaturityIndices[MaturityIndices$Subject == i,])

SubjectNFixed effectsEstimateStd. errorZ valuePr(>|z|)
 180Vocal vs. Postural−1.200.87−1.370.1702
Vocal vs. Locomotor−0.860.78−1.100.2724
Day0.220.063.860.0001
 278Vocal vs. Postural−4.931.54−3.200.0014
Vocal vs. Locomotor−5.351.64−3.260.0011
Day0.360.103.560.0004
 365Vocal vs. Postural−2.151.03−2.090.0362
Vocal vs. Locomotor−4.571.42−3.230.0013
Day0.150.043.550.0004
 478Vocal vs. Postural−3.321.16−2.870.0041
Vocal vs. Locomotor−4.581.38−3.320.0009
Day0.250.073.640.0003
 567Vocal vs. Postural−3.942.01−1.960.0497
Vocal vs. Locomotor−4.511.52−2.970.0030
Day0.210.073.190.0014
 680Vocal vs. Postural−1.010.72−1.400.1607
Vocal vs. Locomotor−1.790.76−2.340.0190
Day0.110.024.42<0.0001
 771Vocal vs. Postural−0.330.82−0.400.6917
Vocal vs. Locomotor0.300.770.390.6989
Day0.100.024.39<0.0001

1.5 Linear regression models to predict contact call duration from locomotor activity. Results are associated with Figure 3B

1.5.1 Population-level linear mixed model results

Formula: Duration ~ LocomotorActivityMean + (LocomotorActivityMean | Day/Subject)

Data: CallParameters

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)0.42230.6498
LocomotorActivityMean1.02541.0126−0.02
Day(Intercept)0.00000.0000
LocomotorActivityMean0.14420.3797NaN
Residual 1.7564 1.3253

Number of obs: 9609, groups: Subject:Day, 188; Day, 38

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)2.291140.05021187.3782545.630<2e-16
LocomotorActivityMean−0.559320.1314936.32386−4.2540.000141

1.5.2 Linear mixed model results per Day

Formula: Duration ~ LocomotorActivityMean + (1|Subject)

Data: CallParameters[CallParameters$Day == i,]

DayNEstimateStd. errorT valuePr(>|t|)
1471−0.460.52−0.880.3794
2369−0.260.43−0.590.5529
33450.010.420.020.9823
4324−0.430.44−0.990.3245
5353−1.340.58−2.330.0203
6403−0.340.39−0.870.3831
7336−0.900.52−1.730.0846
83520.430.510.840.4009
9379−0.510.32−1.570.1164
10258−1.710.31−5.59<0.0001
11297−1.290.31−4.21<0.0001
12336−0.640.33−1.910.0564
13348−1.600.29−5.51<0.0001
14250−0.210.29−0.740.4628
15365−0.500.32−1.540.1243
16113−1.090.52−2.090.0386
17241−0.590.36−1.650.1005
19200−2.130.34−6.25<0.0001
21195−0.760.27−2.820.0053
24217−1.930.28−6.87<0.0001
26255−1.110.29−3.810.0002
28272−0.420.28−1.510.1312
31323−0.520.30−1.730.0851
333120.100.310.310.7555
35222−0.080.33−0.240.8089
38251−0.650.35−1.830.0689
40122−1.210.56−2.170.0318
42269−0.080.25−0.340.7371
45259−0.430.30−1.430.1548
50111−1.010.57−1.770.0793
522810.200.220.880.3786
54177−1.100.35−3.120.0021
561210.810.382.110.0369
592180.330.291.160.2482
611350.730.371.960.0521

1.5.3 Linear mixed model results per Subject

Formula: Duration ~ LocomotorActivityMean + (1 | Day)

Data: CallParameters[CallParameters$Subject == i,]

SubjectNEstimateStd. errorT valuePr(>|t|)
11843−0.830.14−6.06<0.0001
21412−0.900.12−7.57<0.0001
310920.060.150.370.7114
4800−0.310.18−1.740.0817
51688−0.900.21−4.17<0.0001
61253−1.110.26−4.23<0.0001
71521−0.340.20−1.720.0860

1.6 Linear regression models to predict contact call Wiener entropy from locomotor activity. Results are associated with Figure 3B

1.6.1 Population-level linear mixed model results

Formula: WienerEntropyMean ~LocomotorActivityMean + (LocomotorActivityMean | Day/Subject)

Data: CallParameters

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)1.5743881.25475
LocomotorActivityMean1.4534261.20558−0.51
Day(Intercept)1.0164361.00818
LocomotorActivityMean0.0021950.04685−1.00
Residual1.5664481.25158

Number of obs: 9606, groups: Subject:Day, 188; Day, 38

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)−2.50520.191437.1127−13.0851.78e-15
LocomotorActivityMean0.92950.1211130.99847.6753.36e-12

1.6.2 Linear mixed model results per Day

Formula: WienerEntropyMean ~LocomotorActivityMean + (1 | Subject)

Data: CallParameters[CallParameters$Day == i,]

DayNEstimateStd. errorT valuePr(>|t|)
14710.760.312.470.0137
23691.400.245.78<0.0001
33450.680.351.940.0538
43240.760.322.340.0198
53530.600.321.860.0642
64030.400.172.290.0228
73360.820.263.160.0017
83520.900.322.830.0050
93790.550.341.630.1046
102581.180.323.660.0003
112971.170.313.840.0001
123351.400.393.570.0004
133481.530.285.48<0.0001
14250−0.150.34−0.430.6656
153650.660.332.020.0445
161134.621.233.740.0003
172410.620.461.340.1818
192001.680.443.830.0002
211951.720.453.830.0002
242172.940.466.36<0.0001
262551.280.502.590.0102
28271−0.640.59−1.080.2828
31323−0.750.56−1.330.1846
333120.860.631.370.1730
352221.120.591.910.0579
382511.400.721.950.0521
40121−0.551.15−0.480.6334
422690.430.371.160.2465
452591.090.771.420.1562
501111.891.621.170.2464
522810.640.491.300.1951
54177−0.440.47−0.940.3482
561211.050.801.320.1892
592180.670.601.120.2650
611350.290.680.430.6658

1.6.3 Linear mixed model results per Subject

Formula: WienerEntropyMean ~LocomotorActivityMean + (1 | Day)

Data: CallParameters[CallParameters$Subject == i,]

SubjectNEstimateStd. errorT valuePr(>|t|)
118431.250.196.46<0.0001
214121.330.216.48<0.0001
310920.130.140.950.3431
48000.010.150.060.9494
516881.710.179.87<0.0001
612510.830.155.43<0.0001
715200.130.160.820.4142

1.7 Linear regression models to predict locomotor activity during contact calls from contact call type. Results are associated with Figure 3C

1.7.1 Population-level linear mixed model results

Formula: LocomotorActivityMean ~ CallType_Cry0Phee1 + (CallType_Cry0Phee1 | Day/Subject)

Data: CallParameters

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)0.0075840.08709
CallType_Cry0Phee10.0097850.09892−0.98
Day(Intercept)0.0016350.04044
CallType_Cry0Phee10.0026730.0517−0.99
Residual0.0337070.18359

Number of obs: 9609, groups: Subject:Day, 188; Day, 38

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.130090.0120517.8339310.7992.98e-09
CallType_Cry0Phee1−0.071270.0137319.95781−5.1914.47e-05

1.7.2 Linear mixed model results per Day

Formula: LocomotorActivityMean ~ CallType_Cry0Phee1 + (1 | Subject)

Data: CallParameters[CallParameters$Day == i,]

DayNEstimateStd. errorT valuePr(>|t|)
1471−0.080.03−3.000.0029
2369−0.220.04−6.18<0.0001
3345−0.040.03−1.420.1568
4324−0.060.03−2.100.0369
5353−0.050.02−2.160.0327
6403−0.080.03−2.210.0284
7336−0.100.03−3.250.0014
8352−0.060.02−2.510.0132
9379−0.050.02−2.120.0380
10258−0.160.03−4.70<0.0001
11297−0.130.03−4.97<0.0001
12336−0.070.02−3.220.0032
13348−0.170.03−6.28<0.0001
142500.010.020.400.6867
15365−0.100.02−5.03<0.0001
16113−0.130.03−3.700.0003
17241−0.110.03−3.340.0038
19200−0.230.03−6.51<0.0001
21195−0.320.04−8.01<0.0001
24217−0.320.04−7.98<0.0001
26255−0.110.04−2.730.0156
28272−0.020.02−0.980.3512
313230.000.020.120.9076
33312−0.090.04−2.420.0162
35222−0.050.07−0.780.4344
382510.000.03−0.100.9180
40122−0.030.03−0.910.3650
501110.070.100.670.5068

1.7.3 Linear mixed model results per Subject

Formula: LocomotorActivityMean ~ CallType_Cry0Phee1 + (1 | Day)

Data: CallParameters[CallParameters$Subject == i,]

SubjectNEstimateStd. errorT valuePr(>|t|)
11843−0.070.01−7.20<0.0001
21412−0.150.02−9.19<0.0001
31092−0.010.01−0.490.6280
4800−0.020.02−1.000.3172
51688−0.130.01−10.97<0.0001
61253−0.100.01−8.48<0.0001
71521−0.040.01−2.960.0035

1.8 Linear regression models to predict locomotor activity during mature contact calls from age group (early vs late development). Results are associated with Figure 3D

1.8.1 Population-level linear mixed model results

Formula: LocomotorActivityMean ~ AgeGroup_Early0Late1 + (AgeGroup_Early0Late1| Day/Subject)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1,]

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)7.276e-040.0269731
AgeGroup_Early0Late11.774e-030.04211530.39
Day(Intercept)0.000e+000.0000000
AgeGroup_Early0Late11.060e-070.0003255NaN
Residual1.904e-020.1379996

Number of obs: 1846, groups: Subject:Day, 82; Day, 15

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)0.0434920.00616550.1941747.0554.83e-09
AgeGroup_Early0Late10.0424080.0139223.6938563.0470.00561

1.8.2 Linear mixed model results per Subject

Formula: LocomotorActivityMean ~ AgeGroup_Early0Late1 + (1|Day)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1 and CallParameters$Subject == i,])

SubjectNEstimateStd. errorT valuePr(>|t|)
12950.0050.0200.270.7998
2339−0.0210.027−0.790.4474
32740.0160.0240.650.5170
42010.0030.0270.110.9109
52720.1860.0424.410.0009
62620.0460.0251.840.0789
72030.0630.0351.790.0932

1.9 Linear mixed models to predict heart rate percentile during contact calls from contact call type. Results are associated with Figure 4A

1.9.1 Population-level linear mixed model results

Formula: HeartRatePercentileMean ~CallType_Cry0Phee1 + (CallType_Cry0Phee1 | Day/Subject)

Data: CallParameters

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)72.378.507
CallType_Cry0Phee196.649.831−0.53
Day(Intercept)0.000.000
CallType_Cry0Phee110.213.195NaN
Residual611.8824.736

Number of obs: 6215, groups: Subject:Day, 149; Day, 37

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)44.6301.14258.09539.078<2e-16
CallType_Cry0Phee13.3761.45327.9992.3240.0276

1.9.2 Linear mixed model results per Day

Formula: HeartRatePercentileMean ~CallType_Cry0Phee1 + (1 | Subject)

Data: CallParameters[CallParameters$Day == i,]

DayNEstimateStd. errorT valuePr(>|t|)
148−9.376.39−1.470.1495
2252−9.854.22−2.330.0204
312015.524.703.300.0013
42395.423.101.750.0817
51234.544.171.090.2781
62132.953.580.820.4110
71875.655.581.010.3124
82061.864.060.460.6499
9198−2.164.69−0.460.6494
10189−1.674.79−0.350.7280
11148−0.395.45−0.070.9425
12244−9.024.05−2.220.0272
132203.914.390.890.3743
1422614.575.102.860.0094
1528510.223.383.030.0031
171573.457.540.460.6549
19415.789.740.590.5815
2176−1.056.91−0.150.8793
241986.247.250.860.3922
262046.745.251.280.2556
2821336.949.673.820.0003
3123110.275.032.040.0438
33277−2.378.83−0.270.7882
352149.0313.510.670.5044
382251.665.010.330.7423
40112−2.616.57−0.400.6913

1.9.3 Linear mixed model results per Subject

Formula: HeartRatePercentileMean ~CallType_Cry0Phee1 + (1 | Day)

Data: CallParameters[CallParameters$Subject == i,]

SubjectNEstimateStd. errorT valuePr(>|t|)
16771.092.830.390.7003
2454−4.864.23−1.150.2528
310116.422.532.540.0114
47313.403.021.120.2618
51261−0.881.87−0.470.6393
693510.802.554.23<0.0001
711467.682.203.490.0006

1.10 Linear regression models to predict heart rate percentile during mature contact calls from age group (early vs late development). Results are associated with Figure 4B

1.10.1 Population-level linear mixed model results

Formula: HeartRatePercentileMean ~AgeGroup_Early0Late1 + (AgeGroup_Early0Late1 | Day/Subject)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1,]

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)9.717e + 019.857339
AgeGroup_Early0Late13.598e + 0218.967551−0.82
Day(Intercept)2.077e-060.001441
AgeGroup_Early0Late19.461e-060.003076−0.88
Residual5.561e + 0223.580851

Number of obs: 1270, groups: Subject:Day, 58; Day, 15

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)48.132.07428.72223.205<2e-16
AgeGroup_Early0Late16.2773.47836.1961.8050.0795

1.10.2 Linear model results per Subject

Formula: HeartRatePercentileMean ~AgeGroup_Early0Late1 + (1|Day)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1 and CallParameters$Subject == i,])

SubjectNEstimateStd. errorT valuePr(>|t|)
11109.4614.450.650.5524
269−13.2013.26−1.000.4793
3268−14.305.35−2.670.0440
418410.786.821.580.1521
521124.6310.992.240.0521
624510.846.561.650.1424
71832.0310.530.190.8526

1.11 Linear regression models to predict heart rate percentile during locomotor events from age group (early vs late development). Results are associated with Figure 4D

1.11.1 Population-level linear mixed model results

Formula: HeartRatePercentileMean ~AgeGroup_Early0Late1 + (AgeGroup_Early0Late1| Day/Subject)

Data: LocomotorParameters

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)89.149.441
AgeGroup_Early0Late1369.7319.228−0.96
Day(Intercept)61.317.830
AgeGroup_Early0Late161.367.833−1.00
Residual669.6725.878

Number of obs: 1318, groups: Subject:Day, 58; Day, 15

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)44.7173.1618.97414.151.93e-07
AgeGroup_Early0Late18.1164.03920.2512.010.058

1.11.2 Linear mixed model results per Subject

Formula: HeartRatePercentileMean ~AgeGroup_Early0Late1 + (1|Day)

Data: LocomotorParameters[LocomotorParameters$Subject == i,])

SubjectNEstimateStd. errorT valuePr(>|t|)
1106−1.526.40−0.240.8126
238−20.3428.04−0.730.5537
3254−0.508.84−0.060.9569
4187−1.2410.45−0.120.9092
523517.087.802.190.0574
625210.143.512.890.0042
724620.546.463.180.0137

1.12 Linear regression models to predict heart rate percentile from locomotor activity during mature contact calls. Results are associated with Figure 5B

1.12.1 Population-level linear mixed model results

Formula: HeartRatePercentileMean ~LocomotorActivityMean + (LocomotorActivityMean | Day/Subject)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1,]

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)79.088.893
LocomotorActivityMean636.3125.225−0.17
Day(Intercept)11.643.412
LocomotorActivityMean27.745.267−1.00
Residual581.7524.120

Number of obs: 3966, groups: Subject:Day, 144; Day, 37

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)47.5851.09335.73643.538<2e-16
LocomotorActivityMean9.6664.04649.362.3890.0208

1.12.2 Population-level linear mixed model results, data limited to postnatal days 1–30

Formula: HeartRatePercentileMean ~LocomotorActivityMean + (LocomotorActivityMean | Day/Subject)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1 and CallParameters$Day <= 30,]

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)96.57189.8271
LocomotorActivityMean1097.465533.1280−0.64
Day(Intercept)0.11340.3368
LocomotorActivityMean1.18941.09061.00
Residual592.448524.3403

Number of obs: 1705, groups: Subject:Day, 86; Day, 22

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)47.1631.33544.79335.333<2e-16
LocomotorActivityMean16.3176.60331.4462.4710.0191

1.12.3 Population-level linear mixed model results, data limited to postnatal days 31–61

Formula: HeartRatePercentileMean ~LocomotorActivityMean + (LocomotorActivityMean | Day/Subject)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1 and CallParameters$Day > 30,]

Random effects:

GroupsNameVarianceStd.dev.Corr
Subject:Day(Intercept)60.967.808
LocomotorActivityMean421.4020.5280.41
Day(Intercept)29.965.473
LocomotorActivityMean74.548.634−1.00
Residual572.123.919

Number of obs: 2261, groups: Subject:Day, 58; Day, 15

Fixed effects:

EstimateStd. errorDfT valuePr(>|t|)
(Intercept)48.2161.91114.22025.2353.27e-13
LocomotorActivityMean2.6275.24218.8180.5010.622

1.12.4 Linear mixed model results per Day

Formula: HeartRatePercentileMean ~LocomotorActivityMean + (1 | Subject)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1 and CallParameters$Day == i,]

DayNEstimateStd. errorT valuePr(>|t|)
25652.9929.901.770.0820
346−9.9619.25−0.520.6076
47013.7318.410.750.4584
5374.9625.530.190.8471
643−96.7336.48−2.650.0115
7258.0453.550.150.8819
862−24.7919.54−1.270.2096
95462.7422.242.820.0068
106649.7431.191.590.1161
114087.3934.412.540.0153
12115−5.6116.32−0.340.7317
1370−19.2932.42−0.590.5540
1414613.5413.361.010.3125
1514640.1817.612.280.0239
1711018.6622.670.820.4122
1915−23.7126.24−0.900.3827
21599.7625.170.390.6998
241746.5116.620.390.6957
2616839.1812.933.030.0028
2818115.8015.001.050.2938
31159−10.5613.82−0.760.4461
33270−8.9910.22−0.880.3799
35211−3.9811.09−0.360.7202
3817714.3112.221.170.2433
409351.5229.901.720.0882
4226211.598.971.290.1975
4518144.6412.133.680.0003
52241−1.0510.90−0.100.9234
54175−48.9115.80−3.100.0023
5611925.3017.631.430.1540
59121−6.2911.51−0.550.5855
61124−23.3214.14−1.650.1017

1.12.5 Linear mixed model results per Subject

Formula: HeartRatePercentileMean ~LocomotorActivityMean + (1 | Day)

Data: CallParameters[CallParameters$CallType_Cry0Phee1 == 1 and CallParameters$Subject == i,]

SubjectNEstimateStd. errorT valuePr(>|t|)
144721.038.612.440.0150
2351−3.1210.74−0.290.7713
3826−3.565.00−0.710.4770
46070.715.630.130.9002
577531.316.974.49<0.0001
6566−9.679.56−1.010.3120
739418.9310.321.830.0674
https://doi.org/10.7554/eLife.41853.011

Data availability

As part of this full submission, we have uploaded our data to Dryad (https://doi.org/10.5061/dryad.cp75158).

The following data sets were generated
    1. Gustison ML
    2. Borjon JI
    3. Takahashi DY
    4. Ghazanfar AA
    (2019) Dryad Digital Repository
    Data from:Vocal and locomotor coordination develops in association with arousal state.
    https://doi.org/10.5061/dryad.cp75158

References

    1. Adolph KE
    (2008)
    Encyclopedia of Infant and Early Childhood Development
    359–373, Motor and physical development: locomotion, Encyclopedia of Infant and Early Childhood Development, San Diego, CA, Academic Press, 10.1016/B978-012370877-9.00104-3.
  1. Book
    1. Adolph KE
    2. Robinson SR
    (2015)
    Motor development. In: Handbook of child psychology and developmental science
    In: Lerner RM, Liben L, Muller U, editors. Cognitive Processes, 2 (7). New York: Wiley. pp. 114–157.
  2. Book
    1. Blumberg MS
    (2001) The developmental context of thermal homeostasis
    In: Blass E. M, editors. Developmental Psychobiology. Handbook of Behavioral Neurobiology, 13. Boston, MA: Springer. pp. 199–228.
    https://doi.org/10.1007/978-1-4615-1209-7_6
  3. Book
    1. Gesell A
    (1929)
    Infancy and Human Growth
    New York: The Macmillan Company.
  4. Book
    1. Ghazanfar AA
    2. Takahashi DY
    (2017)
    The evo-devo of vocal communication: Insights from marmoset monkeys
    In: Kaas JH, editors. The Evolution of Nervous Systems (2nd). Elsevier Press. pp. 317–324.
    1. Holderied MW
    2. von Helversen O
    (2003) Echolocation range and wingbeat period match in aerial-hawking bats
    Proceedings of the Royal Society of London. Series B: Biological Sciences 270:2293–2299.
    https://doi.org/10.1098/rspb.2003.2487
  5. Book
    1. Hrdy SB
    (2007) Evolutionary context of human development: The cooperative breeding model
    In: Salmon C. A, Shackelford T. K, editors. Family Relationships: An Evolutionary Perspective. Oxford: Oxford University Press. pp. 9–32.
    https://doi.org/10.1093/acprof:oso/9780195320510.003.0003
  6. Book
    1. McGraw M
    (1943)
    The Neuromuscular Maturation of the Human Infant
    New York: Columbia University Press.
  7. Book
    1. Pfaff D
    (2006)
    Brain Arousal and Information Theory
    Cambridge, MA: Harvard University Press.
  8. Book
    1. Rosenblatt JS
    (1976)
    Stages in the early behavioural development of altricial young of selected species of non-primate mammals
    In: Bateson PPG, Hinde RA, editors. Growing Points Ethology. Cambridge: Cambridge University Press. pp. 345–381.
  9. Book
    1. Solomon NG
    2. French JA
    (2007)
    Cooperative Breeding in Mammals
    Cambridge: Cambridge University Press.
  10. Book
    1. Starck JM
    2. Ricklefs RE
    (1998)
    Patterns of development: The altricial-precocial spectrum
    In: Starck JM, Ricklefs RE, editors. Avian Growth and Development. Evolution Within the Altricial Precocial Spectrum. New York: Oxford University Press. pp. 3–30.
    1. Suthers RA
    2. Thomas SP
    3. Suthers BJ
    (1972)
    Respiration, wing-beat and ultrasonic pulse emission in an echo-locating bat
    The Journal of Experimental Biology 56:37–48.
  11. Book
    1. Thelen E
    (1991)
    Motor aspects of emergent speech: A dynamic approach
    In: Krasnegor NA, Rumbaugh DM, Schiefelbusch RL, Studdert-Kennedy M, editors. Biological and Behavioral Determinants of Language Development. Hillsdale, NJ: Lawrence Erlbaum Associates. pp. 339–362.
    1. Tipps ST
    2. Mira MP
    3. Cairns GF
    (1981)
    Concurrent tracking of infant motor and speech development
    Genetic Psychology Monographs 104:303–324.
    1. Vihman MM
    (2014)
    The First Two Years
    Phonological development, The First Two Years, 2nd Edition, Chichester, UK, John Wiley.
    1. Williams H
    (2001)
    Choreography of song, dance and beak movements in the zebra finch (Taeniopygia guttata)
    The Journal of Experimental Biology 204:3497–3506.
    1. Wong JG
    2. Waters DA
    (2001)
    The synchronisation of signal emission with wingbeat during the approach phase in soprano pipistrelles (Pipistrellus pygmaeus)
    The Journal of Experimental Biology 204:575–583.
  12. Book
    1. Yingling JM
    (1981)
    Temporal Features of Infant Speech: A Description of Babbling Patterns Circumscribed by Postural Achievement
    University of Denver.
  13. Book
    1. Zuberbühler K
    (2012)
    Cooperative breeding and the evolution of vocal flexibility
    In: Tallerman M, Gibson K. R, editors. The Oxford Handbook of Language Evolution. Oxford: Oxford University Press. pp. 71–81.

Decision letter

  1. Ofer Tchernichovski
    Reviewing Editor; Hunter College, United States
  2. Ronald L Calabrese
    Senior Editor; Emory University, United States
  3. Franz Goller
    Reviewer; University Utah, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Energetic costs and locomotor constraints on vocal development" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom as a guest Reviewing Editor, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Franz Goller (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Gustison and colleagues recorded and analyzed the development of coordination between vocal and locomotor behaviors in the marmoset monkey. Results show that mature vocalizations develop earlier than locomotion. Further, coordination between vocalization and locomotion develop gradually. Younger infants produce mature-sounding vocalizations only while not moving. Older infants, on the other hand, coordinated their vocalization with movements. Recording heart rate suggest that energetics constraints could potentially explain the late development of coordinated behavior across modalities. However, the reviewers found some substantive concerns, as elaborated below.

Essential revisions:

There are two very major concerns that should be fully addressed before resubmission can be considered:

First, all reviewers concluded (during discussion) that data do not support the major claims authors make about energetic constraints on the development of vocal-locomotor coordination. Therefore, authors should drop out heart rate as a measure of metabolic cost of vocalization. Instead, the revised manuscript should be focused on describing the developmental trajectories with more details, and compare them to relevant literature in other species (which should be presented in Introduction).

Second, authors used sessions (and sometimes even calls) as independent measurements to make statistical inferences. In the revised version, statistics should be redone using each animal as a statistic (i.e., with df = number of monkeys -1). It is essential to test if the phenomenon of different maturation rate can hold using simple paired comparisons. Further, performing the statistics on the fitted spline curves makes it hard to assess if conclusions are sufficiently supported. The dynamic analysis of calls/locomotion coordination should also be supported with similar statistics.

Reviewer #1:

This is an important longitudinal study that compare maturation rate of vocalization to movement and postures. Authors discoveries, if correct, are of general theoretical and practical interest. Of particular interest is the dynamic analysis of movements and calls coordination, showing reversal of patterns over development. Before this manuscript can be published, several concerns about the statistical approach and about conclusions should be addressed:

A) Statistical approach:

1) Using each call as independent statistic is problematic, and (to my knowledge) not an acceptable practice in acoustic analysis. At least for the most critical tests, the n must be the number of monkeys studied, namely 7.

2) Authors should be more careful with their statistical modeling, and take a simple, direct approach whenever possible. In particular, when comparing time courses of calls and movement maturity indexes, author should simply test if those time courses are different. That is, the null hypothesis should be that those are the same. Here all we need is a pairwise comparison within each subject: is the zero crossing the same or not? Such pairwise comparison with n=7 using paired Wilcoxon, or even paired t or just binomial test would be (one can safely reject the null hypothesis (with p=1/2^7) and conclude that the maturation time course of vocalization is indeed faster.

3) The way authors did the test is inappropriate in several ways. First, raising 3 distinct hypotheses (H1= more, H2=same, H3=less) is not necessary, it suffice to look for a significant difference from null. Second, and most critical – where do 124 degrees of freedom comes from? This looks like an obvious pseudo-replication, which resulted with illegitimate power of W = 3936 (!) and an erroneously huge significance (p < 0.0001). Such p values are rarely correct in biology, certainly not with the sample size of the current study.

4) Similarly, the analysis of movement vs. Wiener entropy should be done with each monkey as a statistic, and again, paired comparisons should be done within each animal.

5) The only place where I find it reasonable to use calls as statistics is in the dynamic analysis shown in Figures 4 and 5. In such a fine grain developmental analysis, and after the major effects where established, it makes sense to look at time courses of call movements interactions, considering each event as statistically independent.

B) Conclusions:

This is an important, but still a correlational study – we don't know cause and effect. For example, the statement "infants must overcome the energetic costs of coordinating their vocalizations with body movement" is not fully supported by the data. For example, the marmoset is likely to vocalize while already agitated and this effect might be age related. Further, one might suspect that vocalization might be used to communicate high energy/agitation state.

Reviewer #2:

In this very interesting and well-written manuscript Gustison et al. address the maturation of two motor systems (vocal and locomotor) in marmosets and test a string of clearly defined hypotheses. The authors quantify vocal, postural and locomotory behaviors over complete postnatal development and the resulting datasets are impressive. They use these observations to compute maturity indices for the three behavioral classes, allowing comparison of maturity onsets for the systems separately. The authors convincingly show that young infants only produce adult like calls when not moving and older ones can move while calling. The paper includes several novelties: i) concurrent analyses of both the developmental timescale as well as the real-time (seconds) timescale to investigate energetic tradeoffs during vocalizations. ii) which allows to them study coupling for the first time. The authors show that this coupling exists at the second scale suggesting higher energy investment in calling behavior is required for producing mature calls. The data is well presented with clear figures.

Major comments:

1) The data is used to support the hypothesis that mature contact call production requires an energy investment that increases with age. In the light of energy expenditure during vocalizations this can be explained if the calls get louder. Increasing respiratory effort and thus lung pressure increases predominantly sound amplitude in laryngeally produced sounds. Fundamental frequency only changes slightly with pressure. Wiener entropy as a measure for noisiness is not regularly used outside of the birdsong field, but will increase most likely in chaotic regimes when pressures are very high.

It would be rather straightforward for the authors to quantify the source levels of the vocalizations. The low-frequency cry is not directional but the wave-number can be used to coarsely estimate the directionality of the phee call. (ka = 2pi/λ * a (, where λ is wavelength (v/f) and a is emitter size radius – guessing 2 cm for a marmoset. with k=2pi/(340/7000) and a = 0.02) ka=2.5, which is above 1. Thus the phee call is directional, which makes it more complicated to compute source levels without microphone array. However luckily the microphone was located in the far-field 90 cm above the cage directing down, making vocalization's aimed in the horizontal plane directly comparable.

2) The authors list four empirical data papers in the last paragraph of the discussion that investigated human infant vocal-locomotor coordination. At least Berger et al. and Abney et al. present similar datasets of concurrent posture and vocalization (not at the second timescale). These papers should be acknowledged and placed in perspective up front in the Introduction.

3) The locomotor behavior was basically scored as difference images. Did the authors observe any structural changes in locomotory patterns, such as certain accompanying displays (vocal postures) that matured over development?

Reviewer #3:

General comments:

This manuscript describes a research effort, in which development of locomotor and vocal skills as well as, importantly, their coordination have been studied quantitatively and have been related to energetic aspects. This is a very interesting study that tries to fill a gap in our understanding of the ontogeny of motor system coordination and its energetic background. The presented data partly meet the goals set out in the Introduction (questions 1-3), but especially for questions 2 and 3 fall somewhat short and require clarification.

1) The various sets of hypotheses listed in the results present sets that leave out potential alternatives. Specific suggestions are mentioned below.

2) The assessment of the energetic cost of call production is problematic at multiple levels, from the premises to the estimate and the interpretation. Specifically, the following general points should be considered:

2.1) Several assumptions had to be made to convert heart rate to metabolic activity. We know from critical assessments of the technique that careful calibrations are needed to achieve reliable estimates. That means that the applied conversion factors may not be valid across development within an individual and certainly are problematic when applied across individuals. Whereas this is a general difficulty with the technique that cannot be solved in this study, it should be acknowledged and should lead to more careful conclusions.

2.2) It is not clear to me how baseline levels of heart rate have been determined for call production episodes with and without locomotor activity. Because heart rate changes are very dynamic, this is one of the critical methodological issues, which we face. Are observed changes due to actual changes or differences in baseline levels?

2.3) Rapid changes in heart rate can be caused by motivational and other modulatory factors that have no direct link to energy expenditure of a behavior. This is a major shortcoming of using heart rate to assess the metabolic cost for short-duration behaviors.

2.4) It is not clear to me how the cost of locomotion, which is expected to be several-fold greater than that for vocalization has been disentangled to conclude that production of phee calls imposes a high metabolic cost.

3) The presentation of the results is interwoven by methods (somewhat necessitated by the fact that the method section is at the end) and lots of introductory material and posing of hypotheses. The actual results are then described very briefly, and statistical as well as descriptive presentation is quickly reduced to very derived spline curves. At least to me this generates the feeling that I do not know the results sufficiently to assess whether or not the later conclusions are supported sufficiently.

Specific comments:

Abstract: Is the intended meaning that locomotion uses resources that would otherwise be available for call production? I cannot see that such an energetic bottleneck exists. If a constraint exists, it is more likely that locomotion and call production require either muscle efforts that are difficult to generate simultaneously or some other coordination-related trade-off. The absolute energy level is very unlikely to be the bottle neck here.

Introduction paragraph one: Of course, there are always energetic costs to vocal production. How much energy is required may vary, but as discussed above it is hardly so high that energy levels prohibit vocalization.

Subsection “The vocal system matures before postural and locomotor systems” paragraph two: This description of hypotheses is almost "trivial" and thus does not need to be so elaborate.

Subsection “Mature contact call production and locomotor activity become increasingly coordinated during development”: These hypotheses are based on the same notion that absolute energy levels at the time of vocalization could dictate whether or not calls can be produced. I do not see that such an energy shortage could exist. It is more likely (as alluded to above) that other factors dictate whether or not a call can be produced (e.g., can a particular respiratory muscle generate sufficient force if it also is engaged in locomotion, etc.). Another matter is whether or not an animal engages in a particular behavior because of general energetic condition. These conceptual issues are at the heart of this paper, and, in my opinion, the current reasoning is not realistic in regard to the physiological understanding of metabolic cost and movement.

Subsection “Mature contact call production and locomotor activity become increasingly coordinated during development”: This discussion would read very differently if the constraint were seen more like discussed above.

Final sentence of subsection “Mature contact call production and locomotor activity become increasingly coordinated during development”: This statement may have to be revised if calibrations at different ages were available.

Materials and methods: It is not clear how baseline heart rate was determined. Because metabolic cost of locomotion is not restricted in time to the actual duration of movement, but elevated metabolism continues after the movement, the timing of call production relative to locomotion is very important. This issue is very critical, if the costs of locomotion and vocalization are to be disentangled. It is also not clear how z-scores were calculated to make data comparable over developmental time.

[Editors' note: comments from the first round of re-review follow.]

Thank you for resubmitting your work entitled "Coordination of vocal and locomotor behavior emerges in development as a function of arousal state" for further consideration at eLife. The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below.

In particular:

1) Statistics need improvement. Please provide cross validation by doing shuffle statistics (shuffling subjects across groups). Also, if you choose to continue with the GLM please provide detailed statistical results including within and between subjects results.

2) As reviewer 3 noted, the problem of disentangling metabolism related to locomotion from emotional influences is not fully resolved by using the term 'arousal'. A weaker interpretation of the results could help (in addition to clarifying terminology as suggested).

Reviewer #1:

Authors addressed my comments and the manuscript is much improved but I am still worried about the statistics. Authors now use GLM throughout but they did not specify how dependent repeated observations are treated in the model. judged by the p values, which are all lower than 0.0001, I suspect that p values are based on the dependent call events and movement events as the statistics, rather than the animal. This is a problem. I will need to see the breakdown of within and between subject results, and how p values were computed to confirm, but looking at the raw data, I am almost sure that p values are wrong. Authors should keep in mind that as opposed to R, Matlab is not a statistical software, and it is very easy to get the model wrong. One should use Matlab statistical packages carefully, particularly when extracting p values from complex models, to make sure all assumptions are correct. I still believe it would be much safer to shuffle statistics here: simply shuffle subjects across groups and extract direct p values. Simple is better than sophisticated.

Reviewer #2:

After carefully reviewing the response to reviewers and the new manuscript, I do not have any substantial concerns and recommend the paper for publication.

Reviewer #3:

The manuscript has been improved following the suggestions. The metabolic data and associated argumentation have been taken out. In its place, the heart rate data were used to assess "arousal". "Arousal" is a vague concept and is not defined very well in this manuscript. Whereas I do not see the issues as critical as in the case of metabolism, I am not clear on whether or not the problem of disentangling metabolism related to locomotion from emotional influences has been solved any better here. Perhaps I still do not understand it correctly, but the issue of baseline data has in my opinion also not been resolved satisfactorily. Namely, the interpretation of whether or not and how much heart rate increased depends on the period of baseline measurement. I am not sure I understand what a session is (letter) for which the normalized data are taken. Is it the period before and after as outlined in the manuscript? Certainly, different locomotor activity during these periods will make comparisons very difficult. Because locomotor behavior increases with maturity level, so will the associated heart rates. To what degree it is a change in arousal, can therefore only be said with confidence if the two influences on heart rate can be disentangled.

[Editors' note: comments from the second round of re-review follow.]

Thank you for resubmitting your work entitled "Coordination of vocalization and locomotion emerges in development in association with arousal state" for further consideration at eLife. Your revised article has been evaluated by Ronald Calabrese (Senior Editor), a Reviewing Editor, and a statistician expert. Unfortunately, there are still very serious issues with statistics that preclude publication of your manuscript in eLife at its present form.

Following your disagreement with the reviewers’ criticism about statistics, we sent your manuscript to a statistician expert. The expert concerns are similar to those raised earlier, but are even more substantial. We are worried about the p-values you reported, which seem inconsistent with your sample size and with the raw measures. Therefore, we cannot accept your manuscript, unless we can validate and replicate your results, according to the guidelines provided by the statistician expert.

We would like to clarify that as much as we find your study interesting and potentially valuable, unless the revised version can fully convince the expert that the stats and p-values are correct and easily replicated, we will have no choice but to reject your manuscript.

Reviewer #4:

Major concerns:

1) Authors use LMM which is fine. They say they have infants as random effect and that is fine too, but I wonder how could authors get with 7 monkeys p<0.0001? I looked at the tables with per-animal results. The variability is high in the coefficients per monkey, sometimes with results of opposite directions. Treating them as a sample does not yield the precision near to that claimed.

2) I cannot exclude the possibility that authors may be right, BUT when submitting such an analysis reproducible computing is a must. The statistical model itself should be specified and the entire code made available. I think that giving this information is minimal requirement, but making the code and the data available at this stage so that reproducibility could be checked is already a requirement made by leading journals and here it is crucial.

3) In particular in this case, it is important to know whether the model considered sessions as nested within monkeys or not. It will have great impact, and you cannot tell it from the information authors provided.

4) In the correlation analysis authors use all observations as if they are independent and that is wrong, but since it is used only for creation of scale it does not matter much.

5) Authors also present too many figures only with no data, only summaries: all of 4 and 5B.

6) eLife requires treatment of multiplicity issues; none of that was done.

[Editors' note: comments from the third round of re-review follow.]

Thank you for submitting your article "Vocal and locomotor coordination develops in association with arousal state" for consideration by eLife. Your article has been reviewed by one peer reviewer, and the evaluation has been overseen by a Reviewing Editor and Ronald Calabrese as the Senior Editor.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

We appreciate authors efforts to correct the statistical issues. However, there are still major problems with statistics that authors need to address before the manuscript can be accepted. Please carefully follow the statistician guidance, and in particular correct the models to include individual slopes, and properly adjust p values for multiplicity. Please send us correspondence if statistical significance of the major finding is no longer valid.

Reviewer #4:

Reviewer notes on "Vocal and Locomotor coordination…" by Gustison et al.

I thank the authors for their response to my previous questions. Their effort to make their analysis transparent and reproducible are evident in this version. I hope this version can serve as a model for reproducible research for others.

These compliments do not mean that there are no problems with their analyses,

I shall list three of them:

1) The model they used allows only random intercepts, and not random slopes.

hence the variability from one monkey to the other is not reflected in the calculations of the slope's standard deviation and p-value, which are far too small. In order to enhance replicability of their results, when the experiment will be conducted on another set of monkeys the random slopes inference is important (Then the models should be ~Day + Day/Subject)

The analyses for individuals is appropriate so it lets assess how serious this flaw is.

Sometimes it is very inappropriate: HeartRatePercentileMean~CallType_Cry0Phee1

The estimated same slope is is 3.6 and StanError.96

The individual ones are 1.09, -4.9, 6.4, 3.4, -0.9, 10.8, 7.7 and from these 7 values the standard errors about 2.

Sometimes it will be more reasonable and the conclusion will probably remain unchanged.

The estimated same slope is is.25 and StanError.06

The individual ones are.25,.39,.11,.25,.2,.11,.11, and from these 7 values the standard errors about.04

2) When modelling individual events within the day (from Analysis 1.5 on) all observations within a day (up to 50) are assumed independent.

In such longitudinal data within the day correlations from observation to the next one are expected to be high.

This hampers analyses of individuals as well producing to optimistic standard errors and small p-values.

A simple way out is to construct a daily summary.

3) Adjusting for multiplicity is needed also when trying different models for heart rate, namely its dependency on different variables and in different subsets (1 month 2nd month etc)

With the current p-values not much will change in the conclusion, with new ones I do not know.

https://doi.org/10.7554/eLife.41853.015

Author response

Reviewer #1:

This is an important longitudinal study that compare maturation rate of vocalization to movement and postures. Authors discoveries, if correct, are of general theoretical and practical interest. Of particular interest is the dynamic analysis of movements and calls coordination, showing reversal of patterns over development. Before this manuscript can be published, several concerns about the statistical approach and about conclusions should be addressed:

A) Statistical approach:

1) Using each call as independent statistic is problematic, and (to my knowledge) not an acceptable practice in acoustic analysis. At least for the most critical tests, the n must be the number of monkeys studied, namely 7.

We have revised the analyses to utilize sophisticated techniques - generalized linear mixed effect models - that control for multiple replicating variables (e.g., postnatal day and subject). Our rationale and description of the models are described in responses 2-4.

2) Authors should be more careful with their statistical modeling, and take a simple, direct approach whenever possible. In particular, when comparing time courses of calls and movement maturity indexes, author should simply test if those time courses are different. That is, the null hypothesis should be that those are the same. Here all we need is a pairwise comparison within each subject: is the zero crossing the same or not? Such pairwise comparison with n=7 using paired Wilcoxon, or even paired t or just binomial test would be (one can safely reject the null hypothesis (with p=1/2^7) and conclude that the maturation time course of vocalization is indeed faster.

We understand, appreciate and agree with the reviewer’s point. We wanted to make use of our complex developmental dataset that varies not just by behavioral category, but also by postnatal day and subject. We realize that our original testing regime was inappropriate because it treated each data point as an independent observation. To correct this error in the revised manuscript, we chose to use a sophisticated statistical technique – generalized linear mixed effect models (GLMM) – to make our comparisons. GLMMs are excellent for complex datasets, enabling the user to test for population-level fixed effects while controlling for random effect variables operating on the individual level. Below, we describe the GLMM model that we used and its outcome.

In the GLMM, we include behavior type (vocal vs. postural, and vocal vs. locomotor) as the fixed effect, with postnatal day and subject id as random effects. We used a ‘logistic’ regression function for this GLMM, as it is most appropriate for fitting sigmoid-shaped curves. We had to adapt our maturity index slightly in order to use a logistic GLMM. A logistic regression requires proportion data (i.e., falls within [0,1]). The maturity index in the original manuscript was scaled between -1 and 1. In the revised manuscript, we calculate the maturity index as the proportion of time spent engaged in adult-like behaviors. Thus, the data structure is still the same as before, just shifted to range between 0 and 1.

The outcome of the GLMM is in line with the outcome reported in the original manuscript. The vocal maturation occurs sooner than maturation of postural and locomotor behavior. Because the logistic GLMM analysis requires proportional data, the plots in Figures 2B,D and E have been rescaled. We also updated the significance indicators in plot E.

3) The way authors did the test is inappropriate in several ways. First, raising 3 distinct hypotheses (H1= more, H2=same, H3=less) is not necessary, it suffice to look for a significant difference from null. Second, and most critical – where do 124 degrees of freedom comes from? This looks like an obvious pseudo-replication, which resulted with illegitimate power of W = 3936 (!) and an erroneously huge significance (p < 0.0001). Such p values are rarely correct in biology, certainly not with the sample size of the current study.

We agree with the reviewer that our original analysis technique was fraught with pseudo-replication and therefore inappropriate. We apologize for this. To correct this problem in the revised manuscript, we used a generalized linear mixed effects model. Please see response 2 for a detailed description of our rationale, the model, and changes made to the manuscript.

4) Similarly, the analysis of movement vs. Wiener entropy should be done with each monkey as a statistic, and again, paired comparisons should be done within each animal.

Extending the rationale described above (response 2), we applied generalized linear models to analyses of movement vs. call duration and movement vs. Wiener entropy. We describe these two models in the Materials and methods. Briefly, we include both locomotor activity and postnatal day as fixed effects and infant identity as a random effect. We updated the Results and Figure 3B accordingly. The outcome of the new analyses has the same main outcomes as before -- locomotor activity is negatively associated with call duration and positively associated with Wiener entropy. These models also show that postnatal day is positively associated with call duration and negatively associated with Wiener entropy, which are findings that replicate previous studies.

5) The only place where I find it reasonable to use calls as statistics is in the dynamic analysis shown in Figures 4 and 5. In such a fine grain developmental analysis, and after the major effects where established, it makes sense to look at time courses of call movements interactions, considering each event as statistically independent.

We now use generalized linear mixed models (GLMM) for most analyses. Please see responses 2-4 for details.

B) Conclusions:

This is an important, but still a correlational study – we don't know cause and effect. For example, the statement "infants must overcome the energetic costs of coordinating their vocalizations with body movement" is not fully supported by the data. For example, the marmoset is likely to vocalize while already agitated and this effect might be age related. Further, one might suspect that vocalization might be used to communicate high energy/agitation state.

In the revised manuscript, we use heart rate as a measure of arousal fluctuations rather than attempting to convert it to “energetic cost”. In addition, we removed causal statements when discussing the results. We agree with the reviewer that these statements are misleading given that our study design is observational. For example, we have re-written the Abstract statement mentioned above to “infants gradually improve coordination between vocalizations and body movement through a process that may be facilitated arousal level changes”. We also include a cautionary statement in the Discussion that “However, it is important to note that we are unable to test causal claims directly due to the observational design of this study”.

Reviewer #2:

[…] The paper includes several novelties: i) concurrent analyses of both the developmental timescale as well as the real-time (seconds) timescale to investigate energetic tradeoffs during vocalizations. ii) which allows to them study coupling for the first time. The authors show that this coupling exists at the second scale suggesting higher energy investment in calling behavior is required for producing mature calls. The data is well presented with clear figures.

Major comments:

1) The data is used to support the hypothesis that mature contact call production requires an energy investment that increases with age. In the light of energy expenditure during vocalizations this can be explained if the calls get louder. Increasing respiratory effort and thus lung pressure increases predominantly sound amplitude in laryngeally produced sounds. Fundamental frequency only changes slightly with pressure. Wiener entropy as a measure for noisiness is not regularly used outside of the birdsong field, but will increase most likely in chaotic regimes when pressures are very high.

It would be rather straightforward for the authors to quantify the source levels of the vocalizations. The low-frequency cry is not directional but the wave-number can be used to coarsely estimate the directionality of the phee call. (ka = 2pi/λ * a (, where λ is wavelength (v/f) and a is emitter size radius – guessing 2 cm for a marmoset. with k=2pi/(340/7000) and a = 0.02) ka=2.5, which is above 1. Thus the phee call is directional, which makes it more complicated to compute source levels without microphone array. However luckily the microphone was located in the far-field 90 cm above the cage directing down, making vocalization's aimed in the horizontal plane directly comparable.

We appreciate these insights and we did not consider this approach to measuring sound amplitude as a more robust indicator of respiratory effort. We hope to use it in the future. The current revision, however, adheres to the requirement outlined in the decision letter that was to remove the entire energy expenditure line of argument/analyses. Thus, this comment, while very much appreciated, is no longer applicable to the revised manuscript. If you would still like us to do the analysis, nevertheless, then we’d be happy to give it try. Do let us know.

2) The authors list four empirical data papers in the last paragraph of the discussion that investigated human infant vocal-locomotor coordination. At least Berger et al. and Abney et al. present similar datasets of concurrent posture and vocalization (not at the second timescale). These papers should be acknowledged and placed in perspective up front in the Introduction.

Thank you for these references. We have added these papers to the Introduction in the paragraph where other human-focused studies are now discussed.

3) The locomotor behavior was basically scored as difference images. Did the authors observe any structural changes in locomotory patterns, such as certain accompanying displays (vocal postures) that matured over development?

We did not observe structural changes in the co-occurrence of specific postural-locomotor behaviors with vocalization across development. In an exploratory data analysis, we tested whether phees (mature contact call) were more likely to co-occur with mature forms of posture-locomotion than with immature forms. These analyses were focused on days in which monkeys produced phees and both immature and mature forms of posture (n = 23 sessions) or both immature and mature forms of locomotion (n = 27 sessions). The results suggested that Phee production did not depend on the maturity level of posture (p = 0.70) or locomotion (p = 0.84). Instead, it appeared that phee production in younger infants depended on whether infants were immobile.

Reviewer #3:

General comments:

This manuscript describes a research effort, in which development of locomotor and vocal skills as well as, importantly, their coordination have been studied quantitatively and have been related to energetic aspects. This is a very interesting study that tries to fill a gap in our understanding of the ontogeny of motor system coordination and its energetic background. The presented data partly meet the goals set out in the Introduction (questions 1-3), but especially for questions 2 and 3 fall somewhat short and require clarification.

1) The various sets of hypotheses listed in the results present sets that leave out potential alternatives. Specific suggestions are mentioned below.

To address several of these comments, we chose to drop the “energy expenditure” part of the manuscript. Instead, we focus on heart rate as an estimate of arousal dynamics.

2) The assessment of the energetic cost of call production is problematic at multiple levels, from the premises to the estimate and the interpretation. Specifically, the following general points should be considered:

2.1) Several assumptions had to be made to convert heart rate to metabolic activity. We know from critical assessments of the technique that careful calibrations are needed to achieve reliable estimates. That means that the applied conversion factors may not be valid across development within an individual and certainly are problematic when applied across individuals. Whereas this is a general difficulty with the technique that cannot be solved in this study, it should be acknowledged and should lead to more careful conclusions.

We addressed this comment by using heart rate directly as an indicator of arousal rather than converting it to an estimate of energy expenditure.

2.2) It is not clear to me how baseline levels of heart rate have been determined for call production episodes with and without locomotor activity. Because heart rate changes are very dynamic, this is one of the critical methodological issues, which we face. Are observed changes due to actual changes or differences in baseline levels?

Thank you for this insight. In our study, we normalize heart rate by converting beats/min to percentiles. This normalization method preserves the heart rate changes within a session but allows the dynamic traces to be scaled around the same median (50th percentile) between observation sessions. This is the same method used in a study of arousal dynamics and vocal production in adult marmosets (Borjon et al., 2016). Thus, we aren’t really comparing “baseline” levels of heart rate, we are comparing deviations of heart rate from session “medians”. In the revised manuscript, we describe this better in the Materials and methods and we refer to “median levels” instead of “baseline levels” in the Results.

2.3) Rapid changes in heart rate can be caused by motivational and other modulatory factors that have no direct link to energy expenditure of a behavior. This is a major shortcoming of using heart rate to assess the metabolic cost for short-duration behaviors.

We understand. In the revised manuscript, we used heart rate as an estimate of fluctuating arousal levels rather than of energy expenditure.

2.4) It is not clear to me how the cost of locomotion, which is expected to be several-fold greater than that for vocalization has been disentangled to conclude that production of phee calls imposes a high metabolic cost.

We appreciate the comment. Loud calls such as the marmoset contact call are likely to be very energetically costly for such a small animal (~300 grams in weight) in the same way that it is for other very small animals. Whether it is several fold less than locomotion depends on the duration and mode of locomotion. In any case, for the revised manuscript, the point is moot. We used heart rate as an estimate of fluctuating arousal levels rather than of energy expenditure.

3) The presentation of the results is interwoven by methods (somewhat necessitated by the fact that the method section is at the end) and lots of introductory material and posing of hypotheses. The actual results are then described very briefly, and statistical as well as descriptive presentation is quickly reduced to very derived spline curves. At least to me this generates the feeling that I do not know the results sufficiently to assess whether or not the later conclusions are supported sufficiently.

Thank you for pointing this out. We have removed the original energetic analysis from the revised manuscript, so part of this comment is no longer relevant. In the revised manuscript, we included more mainstream statistical analyses (generalized linear mixed models) so that readers can better understand the results. Some of these models are described above in responses 2-4. In addition to the models described in responses 2-4, we also use GLMMs to supplement the fine-grained dynamic analyses (described in Materials and methods). We believe these models provide more substance to the revised Results sections on the real-time and developmental dynamics between call production, locomotor activity, and arousal levels.

Specific comments:

Abstract: Is the intended meaning that locomotion uses resources that would otherwise be available for call production? I cannot see that such an energetic bottleneck exists. If a constraint exists, it is more likely that locomotion and call production require either muscle efforts that are difficult to generate simultaneously or some other coordination-related trade-off. The absolute energy level is very unlikely to be the bottle neck here.

This is a great point that we did not consider. Indeed, we were focused on the idea that energy was the bottleneck. Sorry to repeat this over and over again but to be complete: In the revised manuscript, we used heart rate as an estimate of fluctuating arousal levels rather than of energy expenditure.

Introduction paragraph one: Of course, there are always energetic costs to vocal production. How much energy is required may vary, but as discussed above it is hardly so high that energy levels prohibit vocalization.

In the revised manuscript, we used heart rate as an estimate of fluctuating arousal levels rather than of energy expenditure.

Subsection “The vocal system matures before postural and locomotor systems” paragraph two: This description of hypotheses is almost "trivial" and thus does not need to be so elaborate.

Thank you for the suggestions. We simplified this section to describe briefly the “null” (overlapping) and “alternative” (non-overlapping) hypotheses. We deleted several redundant phrases (~50 words).

Subsection “Mature contact call production and locomotor activity become increasingly coordinated during development”: These hypotheses are based on the same notion that absolute energy levels at the time of vocalization could dictate whether or not calls can be produced. I do not see that such an energy shortage could exist. It is more likely (as alluded to above) that other factors dictate whether or not a call can be produced (e.g., can a particular respiratory muscle generate sufficient force if it also is engaged in locomotion, etc.). Another matter is whether or not an animal engages in a particular behavior because of general energetic condition. These conceptual issues are at the heart of this paper, and, in my opinion, the current reasoning is not realistic in regard to the physiological understanding of metabolic cost and movement.

Yes, we were wrong (or at least incomplete) in our thinking. We sincerely appreciate you pointing this out. Your insights not only affect the current manuscript (and the elimination of energetics from it) but also how we are thinking about on-going studies in our lab. In the revised manuscript, our interpretations are now focused more on how fluctuating levels of arousal predict motor output and coordination. We refrain from making strong causal statements in the revised manuscript.

Subsection “Mature contact call production and locomotor activity become increasingly coordinated during development”: This discussion would read very differently if the constraint were seen more like discussed above.

This discussion is now focused around using fluctuating arousal levels to predict motor output.

Final sentence of subsection “Mature contact call production and locomotor activity become increasingly coordinated during development”: This statement may have to be revised if calibrations at different ages were available.

Analyses on changes of “energy expenditure” across development have been removed from the revised analysis.

Materials and methods: It is not clear how baseline heart rate was determined. Because metabolic cost of locomotion is not restricted in time to the actual duration of movement, but elevated metabolism continues after the movement, the timing of call production relative to locomotion is very important. This issue is very critical, if the costs of locomotion and vocalization are to be disentangled. It is also not clear how z-scores were calculated to make data comparable over developmental time.

Please see response above for a discussion of what we meant by “baseline” and how we used heart rate percentiles to standardize arousal dynamics across sessions. We agree that dynamic changes in physiological condition are important. For this reason, we continue to include a version of Figure 4, which shows arousal fluctuations surrounding call production and locomotor activity (this analysis is new and described in the Results). Our data suggest that vocal-locomotor behaviors occur during arousal states that have shifted well before (> 10 s) the behaviors occur, and so averages should help to capture these more drastic changes. Given that the manuscript is no longer focused on trying to disentangle the energy costs of behaviors, the nuisances of timing should be less of a concern.

[Editors' note: responses to the first round of re-review follow.]

Reviewer #1:

We respectfully but wholly disagree with this reviewer’s assessment of our statistical approach but we appreciate the opportunity to clarify what we did. Moreover, we’d be happy to subject our analytical approach to the evaluation of anyone the editor views as a more expert statistician than us. Finally, it is worth noting that all our data will be available for anyone to analyze as they wish.

Authors addressed my comments and the manuscript is much improved but I am still worried about the statistics. Authors now use GLM throughout but they did not specify how dependent repeated observations are treated in the model.

We respectfully disagree with the reviewer that we did not specify how dependent repeated observations are treated in the model. As described in the Materials and methods, we include in the regression models the postnatal day and identities of the infants to control for the effect of those variables.

Judged by the p values, which are all lower than 0.0001, I suspect that p values are based on the dependent call events and movement events as the statistics, rather than the animal.

We understand the reviewer’s concern. As the reviewer knows, the p-value just tell us what is the probability that assuming the null hypothesis we would observe a value that is further a part in the distribution; it is arbitrary. Therefore, the reviewer's inference that the small (in his/her opinion) p-value is automatically a consequence of the "wrong" sample size is not correct. For example, if you have a sample of dependent variable that is [1,2,3] and a corresponding sample of a predictors that is [2,3,4] the correlation value is 1 and p-value is smaller than 0.00000001, although the sample size is only 3.

This is a problem. I will need to see the breakdown of within and between subject results, and how p values were computed to confirm, but looking at the raw data, I am almost sure that p values are wrong. Authors should keep in mind that as opposed to R, Matlab is not a statistical software, and it is very easy to get the model wrong. One should use Matlab statistical packages carefully, particularly when extracting p values from complex models, to make sure all assumptions are correct.

To address this concern, we carried out the regression analyses in R software (packages lme4 and LmerTest) instead of Matlab and updated the revised manuscript accordingly. The outcomes of the R models were compatible with the Matlab models—the p-values are correct. Perhaps this reviewer was concerned about the assumptions. To mitigate against this possibility, we repeated the regressions for individual subjects separately and found similar results. Additionally, we provide the data in Dryad so that readers can test the statistics on their own. Finally, we also cite the scholarship related to our statistical approach using logistic generalized linear mixed effect model in R:

Bates, D., Mächler, M., Bolker, B., and Walker, S. (2014). Fitting linear mixed-effects models using lme4. arXivpreprintarXiv:1406.5823.

Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). lmerTest package: tests in linear mixed effects models. Journal of Statistical Software, 82(13).

I still believe it would be much safer to shuffle statistics here: simply shuffle subjects across groups and extract direct p values. Simple is better than sophisticated.

Our analysis is standard in the developmental literature; it is not extraordinarily sophisticated. We disagree that reducing the longitudinal data to a single value by shuffling subjects is a better way. It might look simpler, but one would simply lose information by doing. From our perspective, it does not do justice to the data to reduce a time series to a mean-and-standard deviation.

Reviewer #3:

The manuscript has been improved following the suggestions. The metabolic data and associated argumentation have been taken out. In its place, the heart rate data were used to assess "arousal". "Arousal" is a vague concept and is not defined very well in this manuscript. Whereas I do not see the issues as critical as in the case of metabolism, I am not clear on whether or not the problem of disentangling metabolism related to locomotion from emotional influences has been solved any better here.

We apologize for not explicitly defining “arousal” in the manuscript. It is in fact a very clearly defined concept even to the degree that investigators can study the genetics of arousal (see the work of Donald Pfaff of Rockefeller University). We have used “arousal” as a focus in three of our recent manuscripts without incident (Zhang and Ghazanfar, 2016; Borjon et al., 2016; Liao et al., 2018 PNAS).

In the revised Introduction, we define “arousal” in the following manner: “An animal would be said to exhibit a high arousal state if it is more alert to sensory stimuli, more motorically active and more reactive (Pfaff, 2006). The role of arousal is essentially a means of allocating metabolic energy (i.e., preparing the body for action).”

It is worth noting that “arousal” is separate from “emotions” or “affect”. Both positive and negative affects could be associated with the same level of arousal.

Perhaps I still do not understand it correctly, but the issue of baseline data has in my opinion also not been resolved satisfactorily. Namely, the interpretation of whether or not and how much heart rate increased depends on the period of baseline measurement.

With apologies, we see that our use of the word “baseline” in the revised manuscript is still quite misleading. What we meant in those areas of the manuscript is that the observed real-time and developmental fluctuations in heart rate are outside the 95% threshold of bootstrapped significance tests (i.e., shuffled permutations of the data). In the revised manuscript, we have taken out the word “baseline” completely and now refer specifically to how observed data changes relative to the 95% threshold of the bootstrap significance tests.

I am not sure I understand what a session is (letter) for which the normalized data are taken. Is it the period before and after as outlined in the manuscript?

“Sessions” are 10-minute observations (taken every ~1-3 days) of infants. All behavioral and physiological data presented in this manuscript were collected during these observation sessions and not from before and/or after the observation sessions. In the revised manuscript, we are more careful to refer to sessions specifically as “observation sessions” so that the intended meaning is clearer.

Certainly, different locomotor activity during these periods will make comparisons very difficult. Because locomotor behavior increases with maturity level, so will the associated heart rates. To what degree it is a change in arousal, can therefore only be said with confidence if the two influences on heart rate can be disentangled.

Thank you for this comment. We believe we created some confusion regarding how we re-scaled the heart rate data so that fluctuations can be compared within and across observation sessions. We describe this process in more detail in the Materials and methods section.

Locomotor activity also goes through a re-scaling process within observation sessions, which we describe in more detail in the Materials and methods. It is important to stress here that we are not attempting to compare absolute levels of heart rate and locomotion but rather relative levels of heart rate (percentiles ranging from 0 to 100) and locomotion (smoothed binary ranging from 0 to 1) fluctuations (Appendix 1). In other words, if the heart rate percentile during a contact call is a 100 and the locomotor activity is a 1, that would mean that an infant was mobile with a high heart rate relative to the all data from the same 10-minute observation session.

Another concern this reviewer may have is that any associations between heart rate and locomotor activity during contact calls are a mere by-product of age. To mitigate this possibility, we tested associations between heart rate and locomotion (during phee calls) for one- and two-month old infants separately. We found a positive association for both age groups and report these findings in the revised manuscript:

“Moreover, this positive association characterized infants that were both one month old (1-30 days; n = 1,705 calls, β ± SE = 14.84 ± 4.62, t = 3.21, p = 0.0013) and two months old (31-61 days; n = 2,261 calls, β ± SE = 7.19 ± 3.51, t = 2.05, p = 0.0407). This means that a positive association between locomotor activity and heart rate percentiles was not simply a by-product of age.”

[Editors' note: responses to the second round of re-review follow.]

Reviewer #4:

Major concerns:

1) Authors use LMM which is fine. They say they have infants as random effect and that is fine too, but I wonder how could authors get with 7 monkeys p<0.0001? I looked at the tables with per-animal results. The variability is high in the coefficients per monkey, sometimes with results of opposite directions. Treating them as a sample does not yield the precision near to that claimed.

We understand this reviewer’s concern and we made effort to make the analysis and code more transparent and reproducible. We think that making the data, code, and the model assumptions clear should solve this issue. The reader can now easily check our statistical results and be convinced (or not) by our evidence. The complete answer is in the response to the next comment.

Because this issue with the interpretation of the p-value came out in the previous reviewing round too, and we failed to explain it adequately (partly because we failed to provide the code to the reviewer), we wanted to clarify some points that we think are relevant. Treating each monkey regression coefficient as a sample makes the analysis underpowered. The reason is because the coefficient is not a “single data point”, but the result of a regression line fitted to several data points. As an extreme example, if we have two regression lines and wanted to test whether the coefficients were different or not, we would not say that the sample size is one in each group, because that would preclude any reasonable statistical inference. Obviously, our situation is more complex as we would like to make a generalization to the population. That is why we considered a mixed effect model to increase the statistical power (type II error) of our models. Maybe at more fundamental level, the p-value is not a measure of precision nor of the strength for evidence of the claims. It is simply the probability that we observe a value as extreme as in the data given that the null hypothesis of no effect is correct. The value will depend on the distribution of the alternative hypothesis that is unknown.

2) I cannot exclude the possibility that authors may be right, BUT when submitting such an analysis reproducible computing is a must. The statistical model itself should be specified and the entire code made available. I think that giving this information is minimal requirement, but making the code and the data available at this stage so that reproducibility could be checked is already a requirement made by leading journals and here it is crucial.

We agree with the reviewer that it is crucial that the analyses are more transparent and reproducible. We have addressed this concern in three ways. First, we describe the linear regression models in greater detail in the Data Analysis section. We explicitly state every component of the models (dependent variables, fixed effects, and random effects) and include the model equations used to generate the results in the R script. Second, in Appendix 1, we provide detailed R software output from the population-level regression models (formulas, random effect variance, fixed effect and intercept estimates, standard errors, t values and p values), as well as the output from models run on individual monkeys. Having this appendix will allow interested readers to examine several aspects of the regression models without having to download the datasets and R/Matlab scripts. Third, with the current submission we provide via dryad all data and code needed to run the regression models (R script and.csv data files) and generate the figures (Matlab script and.mat data files).

3) In particular in this case, it is important to know whether the model considered sessions as nested within monkeys or not. It will have great impact, and you cannot tell it from the information authors provided.

We agree with the reviewer that the treatment of random effects as nested or not could have great impact. To avoid any doubt, we followed the recommendation to keep the mixed effect structure maximal (Barr et al., 2013). This means in particular that, whenever possible, we have the variables as both fixed and random effect, keeping the nesting. We re-formulated the majority of regression models to include postnatal day (i.e., sessions) nested within infant subjects (i.e., monkeys) as the random effect. We now include more detailed descriptions of the regression models and their formulas so that it is clear that the nested random effects are used (see response 2). The new regression models had a subtle influence on the statistical outcomes (updated estimates, t-values, and p-values are presented). However, these changes were not large enough to sway the statistical significance at the p < 0.05 level or our interpretation of the models.

Barr DJ, Levy R, Scheepers C, Tilye HJ (2013) Random effects structure for confirmatory hypothesis testing: Keep it maximal. J Mem Lang 68: 255-278. DOI: https://doi.org/10.1016/j.jml.2012.11.001

4) In the correlation analysis authors use all observations as if they are independent and that is wrong, but since it is used only for creation of scale it does not matter much.

We have updated these analyses to use linear mixed models, which mirrors the analysis approach taken in the rest of the manuscript. We describe these regression models in the Data Analysis section, and report the findings in Table 1. The significance outcomes and interpretation did not change by switching over to LMMs from spearman correlations.

5) Authors also present too many figures only with no data, only summaries: all of 4 and 5B.

We adapted Figures 3, 4, and 5 to include data instead of only summaries. In Figures 3B and Figure 5B, we include raw data points in addition to regression slopes and confidence intervals. In Figures 3C, 3D and all of 4, we include additional subplots to illustrate raw data for each individual.

6) eLife requires treatment of multiplicity issues; none of that was done.

In the revised manuscript, we applied a Bonferroni-Holmes correction to adjust p-values in instances where multiple models were used to answer a similar question. Specifically, we applied this correction when examining proportions of time engaged in various behaviors across development (Table 1), associations between locomotor activity (during contact calls) with multiple acoustic parameters, associations between locomotor activity (during contact calls) for different call types, and associations between ANS activity (during contact calls) for different call types. Applying this correction was useful for making p-values more conservative. However, the changes made to p-values were not enough to alter outcomes and interpretation.

[Editors' note: responses to the third round of re-review follow.]

Reviewer #4:

Reviewer notes on "Vocal and Locomotor coordination…" by Gustison et al.

I thank the authors for their response to my previous questions. Their effort to make their analysis transparent and reproducible are evident in this version. I hope this version can serve as a model for reproducible research for others.

These compliments do not mean that there are no problems with their analyses,

I shall list three of them:

1) The model they used allows only random intercepts, and not random slopes.

hence the variability from one monkey to the other is not reflected in the calculations of the slope's standard deviation and p-value, which are far too small. In order to enhance replicability of their results, when the experiment will be conducted on another set of monkeys the random slopes inference is important (Then the models should be ~Day + Day/Subject)

The analyses for individuals is appropriate so it lets assess how serious this flaw is.

Sometimes it is very inappropriate: HeartRatePercentileMean~CallType_Cry0Phee1

The estimated same slope is is 3.6 and StanError.96

The individual ones are 1.09, -4.9, 6.4, 3.4, -0.9, 10.8, 7.7 and from these 7 values the standard errors about 2.

Sometimes it will be more reasonable and the conclusion will probably remain unchanged.

The estimated same slope is is.25 and StanError.06

The individual ones are.25,.39,.11,.25,.2,.11,.11, and from these 7 values the standard errors about.04

Following the reviewer’s recommendation, we changed the generalized and linear mixed model formulas to include random slopes (described in detail in Material and Methods – Data Analysis). These changes to the model formulas resulted in more conservative estimates, standard errors, t-values, and p-values. The main conclusions of our study did not change.

a) The vocal behavior develops faster than posture/locomotion behaviors.

b) Locomotor activity is reduced during mature calls; vocal-locomotor coordination improves across development.

c) Mature calls are produced at elevated arousal; arousal levels are higher during vocal-locomotor coordination.

We are now confident that our statistical output is an accurate representation of the data structure and is fully reproducible.

2) When modelling individual events within the day (from Analysis 1.5 on) all observations within a day (up to 50) are assumed independent.

In such longitudinal data within the day correlations from observation to the next one are expected to be high.

This hampers analyses of individuals as well producing to optimistic standard errors and small p-values.

A simple way out is to construct a daily summary.

Following the reviewer’s advice, we constructed daily summaries and included them in the Appendix. These summary tables include per Day results where appropriate (i.e., for models with Day as a random effect term). The previous version of the manuscript only emphasized Subjects as a random effect term, whereas both Day and Subject are emphasized in the revised manuscript.

3) Adjusting for multiplicity is needed also when trying different models for heart rate, namely its dependency on different variables and in different subsets (1 month 2nd month etc)

With the current p-values not much will change in the conclusion, with new ones I do not know.

Following the reviewer’s recommendation, we’ve adjusted p-values to account for multiplicity issues in models involving heart rate. We still find an association between locomotor activity and arousal level during first month of development.

https://doi.org/10.7554/eLife.41853.016

Article and author information

Author details

  1. Morgan L Gustison

    Princeton Neuroscience Institute, Princeton University, Princeton, United States
    Present address
    Department of Integrative Biology, University of Texas at Austin, Austin, United States
    Contribution
    Data curation, Formal analysis, Investigation, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1162-8966
  2. Jeremy I Borjon

    1. Princeton Neuroscience Institute, Princeton University, Princeton, United States
    2. Department of Psychology, Princeton University, Princeton, United States
    Contribution
    Data curation, Investigation, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9114-3362
  3. Daniel Y Takahashi

    1. Princeton Neuroscience Institute, Princeton University, Princeton, United States
    2. Department of Psychology, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Data curation, Validation, Investigation, Methodology, Writing—review and editing
    For correspondence
    takahashiyd@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4972-001X
  4. Asif A Ghazanfar

    1. Princeton Neuroscience Institute, Princeton University, Princeton, United States
    2. Department of Psychology, Princeton University, Princeton, United States
    3. Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Validation, Writing—original draft, Writing—review and editing
    For correspondence
    asifg@princeton.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1960-7470

Funding

National Institute of Mental Health (T32 MH065214)

  • Morgan L Gustison

National Institute of Neurological Disorders and Stroke (R01 NS054898)

  • Asif A Ghazanfar

National Science Foundation

  • Jeremy I Borjon

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Lauren Kelly for her help with experiments and careful reading and editing of an earlier version of this manuscript. We thank Steve Phelps and Jon Sakata for their feedback on an earlier version of the manuscript. We also thank Yisi Zhang for advice on analyzing heart rates. This work was supported by NIH T32MH065214 (MLG), NSF Graduate Fellowship (JIB), and NIH R01NS054898 (AAG).

Ethics

Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (#1908-18) of Princeton University.

Senior Editor

  1. Ronald L Calabrese, Emory University, United States

Reviewing Editor

  1. Ofer Tchernichovski, Hunter College, United States

Reviewer

  1. Franz Goller, University Utah, United States

Publication history

  1. Received: September 8, 2018
  2. Accepted: July 6, 2019
  3. Accepted Manuscript published: July 16, 2019 (version 1)
  4. Version of Record published: August 6, 2019 (version 2)

Copyright

© 2019, Gustison et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 809
    Page views
  • 117
    Downloads
  • 5
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Immunology and Inflammation
    2. Neuroscience
    Ibrahim T Mughrabi et al.
    Tools and Resources Updated

    Vagus nerve stimulation (VNS) suppresses inflammation and autoimmune diseases in preclinical and clinical studies. The underlying molecular, neurological, and anatomical mechanisms have been well characterized using acute electrophysiological stimulation of the vagus. However, there are several unanswered mechanistic questions about the effects of chronic VNS, which require solving numerous technical challenges for a long-term interface with the vagus in mice. Here, we describe a scalable model for long-term VNS in mice developed and validated in four research laboratories. We observed significant heart rate responses for at least 4 weeks in 60–90% of animals. Device implantation did not impair vagus-mediated reflexes. VNS using this implant significantly suppressed TNF levels in endotoxemia. Histological examination of implanted nerves revealed fibrotic encapsulation without axonal pathology. This model may be useful to study the physiology of the vagus and provides a tool to systematically investigate long-term VNS as therapy for chronic diseases modeled in mice.

    1. Neuroscience
    Shinya Ohara et al.
    Research Article Updated

    The entorhinal cortex, in particular neurons in layer V, allegedly mediate transfer of information from the hippocampus to the neocortex, underlying long-term memory. Recently, this circuit has been shown to comprise a hippocampal output recipient layer Vb and a cortical projecting layer Va. With the use of in vitro electrophysiology in transgenic mice specific for layer Vb, we assessed the presence of the thus necessary connection from layer Vb-to-Va in the functionally distinct medial (MEC) and lateral (LEC) subdivisions; MEC, particularly its dorsal part, processes allocentric spatial information, whereas the corresponding part of LEC processes information representing elements of episodes. Using identical experimental approaches, we show that connections from layer Vb-to-Va neurons are stronger in dorsal LEC compared with dorsal MEC, suggesting different operating principles in these two regions. Although further in vivo experiments are needed, our findings imply a potential difference in how LEC and MEC mediate episodic systems consolidation.