Determination of sample size
Based on restriction of resources and feasibility considerations, we decided to assess a total of N = 60 participants, with 20 participants in each condition (vagus nerve massage, soft shoulder massage, resting control group), prior to conducting the study. To account for the potential exclusion of some participants, we assessed a total of N = 63 participants overall. After the exclusion of n = 3 (reasons see below), we achieved 80% power to detect effects sizes of at least f = 0.41 at a significance level of α = 0.05.
Sixty-three healthy, female participants were recruited at the University of Constance over a three-month period in the summer of 2018. To minimize the risk of negative side effects of the massage (brief discomfort, nausea or dizziness), and effects of potential modulators of the PNS, exclusion criteria involved past, or current orthopaedic problems (e.g. scoliosis, discopathy, spinal fracture), whiplash, thyroid or cardiac disease, intake of blood-thinning or anticoagulant medication, smoking and left-handedness. These criteria were emphasized on the recruitment flyers, and participants confirmed the compliance to these criteria in written form at the beginning of the experimental procedure. Interested male participants were not invited to control for possible sex effects in the massage giver/receiver interaction, since experimenters were all female. Eligible participants were assigned to either a vagus nerve massage (VNM), a soft shoulder massage (SSM), or a resting condition (resting control group, RCG) by the experimenters using block randomization (block size of three) in a parallel design to assure equivalent sample sizes among groups. Using a block size of three, it was determined prior to the recruitment of participants which experimental condition was tested on which day and time (no double-blind approach). Subsequently, participants signed up blindly to a session and were not informed about their allocation to a certain condition prior to the debriefing. Participants were asked to withhold from consuming any caffeine-containing beverages 2 h prior to testing. Three participants had to be excluded due to technical problems (n = 1), lack of German skills (n = 1), or disruption of the experimental procedure (n = 1). A flow diagram visualizing the sample is depicted in Fig. 3.
The final sample of 60 participants (meanage = 22.76 years, SDage = 3.43, rangeage = [19, 36] years) consisted of students (n = 57; 95.00%), employed (n = 1; 1.67%), self-employed (n = 1; 1.67%), and unemployed healthy adults (n = 1; 1.67%). 61.67% (n = 37) were currently in a relationship (mean duration 2.47 years). The majority (n = 54 of n = 59) stated that they identified themselves as (predominantly) heterosexual (n = 1 participant did not answer that question).
Experimental sessions were scheduled to start between 7:00 a.m. and 7:00 p.m. and lasted for approximately 50 min. The sampling took place in the same room throughout the entire study. The detailed study procedure is depicted in Fig. 4.
During an introduction, participants gave written informed consent, confirmed the compliance to the exclusion criteria in written form, and put on the chest strap for continuous heart rate (HR) monitoring. After that, participants sat down for an acclimatization period, filled in questionnaires and sat quietly for 5 min to obtain a physiological baseline. They were instructed to sit upright with both feet flat on the floor, and knees at a 90° angle. This was followed by the completion of a hand power test and a fine motor skills task (5 min) which are not analysed in the course of this work. Directly afterwards, a standardized 10-min intervention of either a VNM, a SSM, or a resting condition was carried out during which participants remained seated (RCG). After the completion of a creative task, and a post assessment of the hand power test, and the fine motor skills task (5 min) which are not analysed in the course of this work, participants completed the Social Touch Questionnaire at their own pace. This questionnaire period was considered as a recovery period.
Participants were asked to rate their current mood on the Affect Grid (described in detail below) at four time points during the study procedure and their perceived stress level on a modified version of the Perceived Stress Questionnaire at two time points. At the end, participants were thanked, debriefed and received compensation for participation.
Experimental and control conditions
During the 10-min intervention, all participants were asked to sit comfortably, place their head on a massage cushion on the table, and close their eyes. Participants were informed that the massage would be interrupted immediately, as soon as they indicated any negative short-term consequences (brief discomfort, nausea or dizziness). The head was aligned as straight as possible to avoid a one-sided stimulation. Whenever questions occurred, conversations were restricted to a minimum; else the intervention phase was carried out in silence. The VNM group received a standardized head/neck massage that was applied with moderate pressure. The stimulation included stroking and twisting the trapezius muscle by gripping between the trapezius and sternocleidomastoid muscle, below which the main strand of the vagus nerve runs, to reach deeper regions. Further, the VNM involved stroking the muscles below the base of the skull from the spine to the back of the ears. In contrast, the SSM was applied by stroking and softly touching the neck and shoulder area. The RCG did not receive any physical contact. The SPs including detailed instructions for the VNM and the SSM can be obtained from the Supplementary Information.
High frequency heart rate variability (HF-HRV)
RR intervals were collected continuously at a sampling rate of 1,000 Hz using a Bluetooth low energy Polar H7 heart rate sensor on a two-electrode chest strap (Polar, Finland) and the application Heart Rate Variability Logger for iOS34. RR intervals differing more than 20% from the previous interval were discarded by the application to correct for artefacts caused by noise or motion. HF-HRV of the events of interest (baseline, intervention and recovery) was calculated using the R package RHRV35.
HF-HRV reflects the temporal variation between successive heart beats that occurs in the respiratory frequency band31,36,37. It is one of the most commonly used vagally-mediated heart rate variability parameters38 which is generally accepted to validly reflect parasympathetic activity31,36. This has been shown in a number of foundational studies39 showing that scopolamine, a pharmacological block of the PNS, completely eliminates HF-HRV. In contrast, beta-adrenergic blockers, which block the sympathetic nervous system, have no influence on HF-HRV. Thus, a SP triggering a physiological relaxation response could be identified using HF-HRV as a biomarker of PNS activity.
For the extraction of HF-HRV, we employed an R script (in-house) that automatically derives 5 min of each subject’s baseline, 10 min of each subject’s intervention, and 2.5 min of each subject’s recovery phase. These segments were split up into 2.5-min intervals (with exception of the recovery phase segment) and the instantaneous heart rate signal was extracted for each segment. Artefacts were detected automatically by a filtering algorithm provided by RHRV that uses adaptive thresholding to reject beats whose RR value differ from preceding and subsequent beats more than a threshold value. The algorithm further removed values that are not within an acceptable physiological range35. In addition to that, RR intervals were visually screened, and ectopic beats were removed. After removing artefacts, the heart rate signal was interpolated at a sampling frequency of 4 Hz. In the subsequent frequency analysis, mean HF-HRV over the 2.5-min periods was calculated using 60 s segments with a shift of 30 s, and a fixed frequency bandwidth of 0.15–0.4 Hz. Since the duration of the recovery time (the time during which participants completed the Social Touch Questionnaire) was not standardized temporarily across participants, and to allow for a better comparability between subjects, we used a fixed window (2.5 min) for HF-HRV calculations. Since HF-HRV calculations are considered unreliable when the time period of RR assessment is too short, participants with a total recovery time interval shorter than 2.5 min (n = 3 of VNM, n = 1 of SSM, n = 1 of RCG) were assigned no HF-HRV value for this period. Finally, to account for the violation of the normality assumption, we calculated the natural logarithm of the resulting HF-HRV values and used them as an index of HF-HRV in all following statistical analyses. For readability reasons, HF-HRV is used to refer to ln(HF-HRV) throughout the manuscript.
Beside frequency domain heart rate variability components, time domain parameters such as the root mean square of successive differences (2) have been used frequently to describe sympatho-vagal balance of the autonomic nervous system31. In the course of the processing described above, RMSSD was calculated analogously to HF-HRV. Since we were particularly interested in changes in PNS activity, we will focus on the effects of HF-HRV in the following. RMSSD data can be obtained online at the Open Science Framework project associated with this work.
Subjective levels of relaxation
Psychological relaxation was measured using the Affect Grid40, which assesses affective state along the two dimensions displeasure/pleasure, and arousal/sleepiness. The instrument has adequate reliability and validity40,41, and has been used in in different experimental settings, such as longitudinal ecological momentary assessment42, or experimental designs43,44,45. The score for each dimension ranges from 1 to 9. To analyse subjective levels of relaxation, we multiplied the pleasure and arousal scores to receive a single-item score, with higher scores indicating higher levels of subjective relaxation (range 1–81).
Subjective levels of stress
Subjective stress levels were assessed using the revised, 20 item German version of the Perceived Stress Questionnaire (PSQ46), a widely used instrument in research on stress and well-being47. The questionnaire shows satisfactory reliability46, and validity48. We modified the temporal frame by asking participants to rate how strongly an item applied to them on a 4-point Likert scale (1 = does not apply at all, 4 = does entirely apply) in current moment. An overall sum score was calculated after inversing the items of the scale “joy”. The higher the score, the more stress subjects perceived within the situation (range 20–80).
Social Touch Questionnaire
To assess the general attitude towards touch, participants completed the Social Touch Questionnaire (STQ49), a 20 item scale asking participants to rate how strongly an item applies to them on a 5-point Likert scale (0 = not at all, 4 = extremely). The instrument is frequently used in research on social touch and related areas (e.g.50) and is considered a reliable and valid instrument49,51. A sum score was calculated (range 0–80), with lower scores indicating a more positive attitude towards social touch.
Statistical analyses were conducted in R version 3.5.352 using RStudio version 1.1.46353 with the packages ez54, nlme55, and car56, and in Jeffreys’s Amazing Statistics Program57 version 0.11.1. Graphs were created using the R package ggplot258. If required, normality assumption was checked using the Shapiro–Wilk Normality test and histograms of model residuals, and homogeneity of variance assumption was checked using Levene’s test. The level of significance was set at α = 0.05. Overall, n = 19 of the VNM, n = 22 of the SSM, and n = 19 of the RCG were included in the following analyses.
One-way ANOVAs (type II) with Experimental Condition (VNM, SSM, RCG) as independent variable were performed to ensure equal distribution of age, heart rate baseline (mean of the first and second part of the heart rate baseline measurement), and STQ across the three groups. Furthermore, Pearson’s Chi-squared tests were conducted to test whether the variables sexual orientation [(predominantly) heterosexual/other), and session start (morning/afternoon)] were equally distributed across the experimental groups. Variables that were not equally distributed across the groups were considered as potential confounds, and were controlled for in subsequent analyses.
First, we performed mixed analyses of variance (ANOVAs) with Experimental Condition (VNM, SSM, RCG) as between-subject factor, Time as within-subject factor and the outcomes subjective stress and subjective relaxation to assess whether the different interventions induced different changes in subjective stress, and relaxation levels. Greenhouse–Geisser corrections were applied whenever the sphericity assumption was violated according to the Maulchy’s test for Sphericity. Partial eta squared (partial η2) was used as indicator of effect size59,60. In addition, Bayes factors (BF10) were calculated in comparison to the null model, indicating the likelihood of the data given the alternative hypothesis is true divided by the likelihood of the data given the null hypothesis is true61,62.
Because of five missing values in HF-HRV during the recovery period and high inter-individual variations in HRV31, a growth curve approach was performed to assess changes in HF-HRV over time, estimating within-subject trajectories of HF-HRV over time on level 1, individual differences on level 2, and differences caused by experimental condition on level 3. To test the hypothesis that receiving a massage (VNM and SSM condition) leads to greater increases in HF-HRV compared to rest (RCG), experimental condition was entered as a dummy variable called Touch Effect (RCG = 0, SSM = 1, VNM = 1). In the baseline model, HF-HRV was predicted from the intercept. After that, random intercepts across individuals (random intercept model) were incorporated. Next, a linear fixed effect of Time (Time model) was included. Thereafter, random effects of Time over people (random slope model) were incorporated. Then, a quadratic trend of Time (Time2 model) and a cubic trend of Time (Time3 model) were added successively as orthogonal predictors, as HF-HRV was expected to show a curvilinear change over the course of the experiment. After that, a first-order autoregressive covariance structure was added to the model (random intercept, fixed slope model with covariance structure). After adding the condition effect Touch Effect (condition model), the interaction terms Time × Touch Effect, Time2 × Touch Effect, and Time3 × Touch Effect were entered into the model (interaction model). Resulting changes in overall model fit by means of the log-likelihood value were compared using an ANOVA and the final model was evaluated. As an estimate of effect size in nested models, we calculated the likelihood ratio of the condition model and the model with random intercept, fixed slope, and covariance structure, and the likelihood ratio of the interaction model and the model with random intercept, fixed slope, and covariance structure. As such, the likelihood ratio represents how much more likely the data is under the given model compared to the less complex model.
In order to calculate standardized regression coefficients (which are interpretable as effect sizes), we z-standardized the HF-HRV values before conducting the analyses.
To test whether the VNM leads to a higher HRV increase compared to the SSM, we excluded the RCG, and again entered experimental condition as a dummy variable called Vagus Effect (SSM = 0, VNM = 1) and subsequently ran the same analysis as described above.
Ethics approval, written informed consent, and compensation
The study was approved by the Ethics Committee of the University of Constance, Germany, prior to its conductance (IRB statement 12/2017), and was carried out in accordance with the ethical standards of the Declaration of Helsinki. All participants gave written informed consent prior to participation and received financial compensation (€10). The study is officially registered at the German Clinical Trials Register (https://www.drks.de, trial number DRKS00021162, date of registration 27/03/2020).