Effects of a short residential thermal spa program to prevent work-related stress/burnout on stress biomarkers: the ThermStress proof of concept study.

Introduction

Stress at work is both a public health issue and an economic issue.[1,2] The main consequences of stress on health are increased cardiovascular morbidity and mortality.[3] Stress particularly affects autonomic nervous system functioning.[4] One of the main cardiovascular health effects of stress is deregulation of the sympathovagal balance.[5] The most common measures of autonomic nervous system activity are heart rate variability (HRV) and electrodermal activity.[6] Conveniently, there are non-intrusive and pain-free measures of HRV and electrodermal activity.[7] Some genetic polymorphisms linked with stress, such as polymorphism of angiotensin converting enzymes[8-11] or of serotonin,[12-14] are also associated with cardiovascular risk. Stress also causes arterial ischemic pathology[15] via complex mechanisms involving changes in arterial endothelial and microvascular atherosclerosis.[15] These microvascular changes are linked to systemic inflammation caused by stress,[16] and changes within the hypothalamic–pituitary axis (e.g. dehydroepiandrosterone sulphate (DHEAS),[17,18] cortisol,[17,18]) and at a central level (e.g. neuropeptide Y,[19] brain-derived natriuretic factor (BDNF),[20]). Stress can also lead to obesity through inappropriate eating behaviours.[21] Stressed individuals also find it difficult to lose weight.[22] Possible biological mechanisms that link obesity and stress involve the action of stress on leptin,[23] an anorectic hormone secreted by adipose tissue and proposed as a biomarker of stress.[24] Stress, particularly stress at work, is responsible for multiple somatic diseases.[25,26] For example, long-term mental stress can have negative effects on bone tissue, leading to osteoporosis and increased bone fracture risk.[27] Therefore, an assessment of the association between bone parameters and stress is needed. Moreover, our program may modify bone parameters and decrease bone fracture risk.[6,28] More interestingly, peripheral quantitative computed tomography (pQCT) can be used to measure muscles and intra- and inter-muscular fat content,[29-31] which can also be modified by our program[32] and may be biomarkers of stress. Many stress biomarkers are secreted by adipose tissue.[33] Moreover, it has been demonstrated that limb composition reflects total body composition, particularly fat.[29-31]

The psychological consequences of stress are numerous and include dissatisfaction, anxiety, depression and burnout.[26] Non-pharmacological methods of stress management include psychological interventions, physical activity[34,35] and therapy for stress-induced eating disorders.[36] Several psychological interventions have demonstrated positive effects on work-related stress, such as cognitive-behavioural therapy,[37] acceptance and commitment therapy[38,39] and mindfulness.[40,41] The benefits of any type of physical activity on the physical and mental health of individuals at any age are indisputable.[42]

A recent systematic review highlighted the effects of balneotherapy (bathing in mineral springs) and spa therapy on stress biomarkers. This review of 15 studies (684 subjects) emphasised the potential benefits of such treatment on cortisol levels.[43] In France, there are five spa resorts that specialise in the treatment of psychosomatic disorders. The usual duration for a residential thermal program in these resorts is 3 weeks. Previous research has demonstrated the positive effects of a 3-week spa therapy intervention on burnout.[44] However, a shorter residential thermal spa program may be more compatible with a professional context because of the availability of individuals. In addition, it would be useful to focus on work-related stress prevention (i.e. before the state of burnout).

Objective

This ThermStress proof of concept study was developed to provide a new approach to the management of stress disorders. The main aim of the study is to assess the ability of a short residential spa program to increase HRV and manage work-related stress.

Methods

This 1-year study with repeated measures on six occasions (inclusion (one to five months prior to the beginning of the spa program – M-1 to M-5), 6 days before the start (D-6), at the start (D0), at the end of the spa program (D6), at 6 months (M6) and at 12 months (M12)), will allow us to explore the effect of a specific short-term residential thermal spa program on the prevention of occupational burnout and work-related stress in workers, focusing on measures of well-being and cardiovascular morbidity. Each participant will undergo a 6-day residential spa program comprising psychological intervention, physical activity, thermal spa treatment, health education and eating disorder therapy and a follow-up. Each participant will be his/her own control and we will compare variations in parameters between intervention/control periods (short-term variations: D6/D0 and D0/D6; long-term variations: inclusion/D0 compared with D0/M6–M12) (Figure 1).

Figure 1.

The ThermStress protocol. M-5: 5 months before the study; M-1: 1 month before the study.

HRV: heart rate variability.

Selection criteria

Volunteers will be workers aged 18 years to retirement, with a stress visual analogue scale score greater than 50 mm, who wish to follow a residential spa thermal program to manage work-related stress. Participants must have had a stable weight during the last 3 months and have no uncontrolled cardiac, hepatic, renal or endocrine diseases.[32] In compliance with human ethics guidelines, participants must be covered by social health insurance and must sign consent forms.

Potential participants will be excluded if major treatment and/or protocol deviations are observed. Other exclusion criteria are drugs and medical conditions that significantly affect the primary outcome of HRV (e.g. alpha or beta-blockers; arrhythmias or conduction disorders such as bundle branch block, atrioventricular heart block).[45]

The rationale for the sample size estimation is based on HRV, which is a biomarker of both stress and morbidity/mortality.[46,47] In particular, within multiple parameters of HRV, we considered log low frequency (LF)/high frequency (HF) HRV for sample size calculation because it is the parameter that traditionally represents sympathovagal balance, together with global HRV (see description of LF/HF below in the description of the primary outcome section).[4,46] A log LF/HF with low values is associated with good autonomic nervous system adaptation. Based on previous work,[48] we assumed a statistical power greater than 80%, a two-sided type I error of 5% and an individual correlation coefficient set at 0.5 (each subject as his/her own control). Based on the results of a pilot study, we assumed that our response will promote an absolute decrease of the log LF/HF by 20% (with a standard deviation of 45%).[49] Therefore, we need to include n = 43 subjects. Finally, to take into account lost to follow-up, it is proposed to include 56 subjects.

Participants will be male and female volunteers aged 18 years to retirement. Following ethics committee approval, and based on our calculation, 56 volunteers will be enrolled to account for potential dropouts. All participants will be given written information regarding the project and will have to sign consent forms before enrolment. Workers will be recruited by the spa centres using flyers and website announcements. Participants will also be recruited using advertisements provided by the CHU of Clermont-Ferrand or by the Clermont Auvergne University.

Participants will attend a short-term residential spa program lasting 6 days, which combines psychological intervention, physical activity, thermal spa treatment, health education and eating disorder therapy. Physical activity will be diverse (endurance, strength, circuit training) to offer a personalised program compatible with the desires of each participant. The benefits of all types of physical activity on the physical and mental health of individuals of any age are indisputable. Psychological interventions will include various validated therapies for work-related stress, such as cognitive-behavioural therapy, acceptance and commitment therapy and mindfulness. Details on intervention components are shown in Table 1 and the agenda of the intervention program is shown in Table 2.

Table 1.

Description of intervention components.

Intervention	Dose	Description
Spa treatment	6 × 120 min	Spa treatment is a combination of spa bath (i.e. thermal mud, jets), spa shower (i.e. hydromassage) and body massage
Psychological workshop	6 × 120 min	Cognitive-behavioural therapy, acceptance and commitment therapy, and mindfulness
Nutrition workshop	4 × 90 min	Participants must understand: concept of energetic balance, calculation of energy intake, calculation of energy expenditure (of physical activity, thermic effect of food), effects of stress on energy intake (quantity and type of intake), effects of physical activity on food intake regulation, cooking sessions (to cook healthy meals)
	+1 breakfast	During the breakfast session, participants will be provided with information and advice related to the preparation of a healthy breakfast
Supervised physical activity	5 × 120 min	Supervised physical activity will be delivered every day by a qualified exercise physiologist. Physical activity will be delivered as circuit training and will include resistance and aerobic training
Biomarkers of Stress workshop	1 × 90 min	Understanding the purpose of biological changes owing to stress and burnout

Table 2.

Agenda of the intervention program.

	Sunday	Monday	Tuesday	Wednesday	Thursday	Friday	Saturday
7h30–8h00		Breakfast	Breakfast	Breakfast with dietician	Breakfast	Breakfast	Breakfast
8h15–10h15		Spa treatment	Spa treatment	Spa treatment	Spa treatment	Spa treatment	Spa treatment
10h30–12h30		Psychological workshop	Psychological workshop	Psychological workshop	Psychological workshop	Psychological workshop	Psychological workshop
12h45–13h45		Lunch	Lunch	Lunch	Lunch	Lunch	Lunch
13h55–15h25	*	Nutrition workshop	Nutrition workshop	Biomarkers of Stress workshop	Nutrition workshop	Nutrition workshop
15h25–15h30	*	Collation[†]	Collation[†]	Collation[†]	Collation[†]	Collation[†]
15h30–17h30	*	Supervised physical activity	Supervised physical activity	Supervised physical activity	Supervised physical activity	Supervised physical activity
17h30–19h30	Free time	Free time	Free time	Free time	Free time	Free time
19h30–20h30	Dinner	Dinner	Dinner	Dinner	Dinner	Dinner

Arrival of participants at their desired time on Sunday afternoon.

Collation is any kind of healthy combination, such as fruits/yogurt or any other low glycaemic index carbohydrates and protein sources.

After the intervention phase of the study, participants will undergo a 1-year at-home follow-up.

Each participant will undergo a battery of tests (described below). Data collection will be performed six times, as previously described, except for dual-energy X-ray absorptiometry (DXA) and pQCT, which will be performed at inclusion and after 12 months, and cardiac remodelling and function, which will be performed at inclusion and after 6 months. Table 3 outlines all the study outcomes.

Table 3.

Outcomes.

Variables	Type of measure	Measurement modality	References
Biomarkers of stress and cardiovascular risk	Skin conductance	Wrist band electrodes – Movisens	[7]
	Blood flow velocity	Laser speckle contrast imaging	[78,79]
	Myocardial longitudinal strain	Speckle-tracking echocardiography	[85,86]
Genetic polymorphisms	Polymorphism of angiotensin converting enzymes	Blood cells	[8–11]
	Polymorphism of serotonin	Blood cells	[12–14]
Demographics	Age, gender, qualification, personal work status, ethnicity, life and occupational events	Questionnaire	[4]
Clinical measurements		Height, weight, blood pressure, heart rate, waist circumference
Body composition	Muscle mass, fat mass, bone structure	Impedance meter	[80]
		Densitometry X-ray absorption	[6]
		Peripheral quantitative computed tomography	[73–75]
		Quantitative ultrasound	[81]
Psychology and quality of life	Depression	Hamilton Anxiety Rating Scale	[62]
	Anxiety	Hamilton Anxiety Rating Scale	[62]
		State–Trait Anxiety Inventory	[60,61]
	General health	General Health Questionnaire (SF-12)	[71]
	Stress, fatigue, sleep	Visual analogue scale	[4]
	Burnout	Maslach Burnout Inventory	[59]
	Mindfulness	Freiburg Mindfulness Inventory	[22,23]
	Coping	Brief COPE Questionnaire	[63]
	Emotions	Emotion Regulation Questionnaire	[64]
	Perception of work	Karasek’s Job Content Questionnaire	[65]
	Self-efficacy	Perceived Self-efficacy Scale	[66]
	Addiction to work	Work Addiction Risk test	[67]
	Alexithymia	Toronto Alexithymia Scale	[68,82,83]
	Illness perception	Brief Illness Perception Questionnaire, adapted for stress at work	[69,84]
	Metacognition	Metacognitions Questionnaire-30	[70]
Lifestyle	E.g., smoking, alcohol consumption, coffee consumption, food intake	Questionnaires	[4]
	Physical activity	Recent Physical Activity Questionnaire	[72]
Alloplastic load	HbA1c, HDLc and LDLc, TG	Ethylenediaminetetraacetic acid tube	[6]
	Cortisol	Dry tube, serum isolation and deep-freezing	[17,85]
	DHEAS	Dry tube, serum isolation and deep-freezing	[17,85]
	Leptin	Dry tube, serum isolation and deep-freezing	[6]
	BDNF	Dry tube, serum isolation and deep-freezing	[20,86,87]
	Pro-inflammatory cytokines: IL-1β, IL-6, IL-1, TNFα	Dry tube, serum isolation and deep-freezing	[6]
	NPY	Dry tube, serum isolation and deep-freezing	[19]
	Telomere length	Blood; Southern blot or polymerase chain reaction analyses	[88]
	Heart rate variability	Holter	[89]
	Skin conductance – electrodermal activity	Wristband electrodes	[89]

A total of 25 mL of blood will be sampled (i.e. 100 mL in 1 year). HbA1c: haemoglobin A1c; HDLc: high-density lipoprotein cholesterol; LDLc: low-density lipoprotein cholesterol; TG: triglycerides; DHEAS: dehydroepiandrosterone sulphate; BDNF: brain-derived neurotrophic factor; IL: interleukin; TNFα: tumour necrosis factor α; NPY: neuropeptide Y.

To assess the ability of a short residential spa program to manage work-related stress and increase HRV, a biomarker of both stress and morbidity/mortality.46,47

HRV parameters will be assessed over 26 hours using a heart rate transmitter belt simply positioned on the chest, with a 26-hour recording time, a beat per minute within a 25 to 240 range and respiratory rate within a 3 to 70 range (Zephyr™ BioHarness™ BT, Zephyr Technology, Annapolis, USA). The HRV data will be examined according to the recommendations of the European Society of Cardiology and the North American Society (Task Force).[50] HRV will be explored in time and frequency domains.[51] A methodology developed by our team will also be used.[52] Premature atrial and ventricular beats will be automatically discarded and visually checked. In the time domain, R-R intervals, standard deviation of R-R intervals, square root of the mean squared difference of successive R-R intervals (rMSSD), and number of adjacent N-N differing by more than 50 ms divided by the total number of N-N intervals (pNN50) will be analysed. The rMSSD and pNN50 are associated with HF power and hence parasympathetic activity. In the spectral domain, we will analyse LF power (0.04–0.15 Hz), an index of both sympathetic and parasympathetic activity, and HF power (0.15–0.4 Hz), representing the most efferent vagal (parasympathetic) activity to the sinus node. Very low frequency (VLF; 0.003–0.04 Hz) partially reflects thermoregulatory mechanisms, fluctuation in activity of the renin–angiotensin system and the function of peripheral chemoreceptors. LF and HF will also be assessed in normalised units (nu), which represent the relative value of each power component in proportion to the total power minus the VLF component. Thus, LFnu and HFnu are considered to represent best sympathetic and parasympathetic activity, respectively.[53] The LF/HF ratio (i.e. the sympathovagal balance) will also be calculated.[46]

Secondary outcomes will be (1) to demonstrate an improvement in perceived stress and other stress biomarkers following the short residential spa program, (2) to study the influence of genetic polymorphisms on stress and on the response to our stress management program, (3) to examine the relationship between subjective (questionnaire responses) and objective (biomarker) variables, (4) to propose a salient biomarker or a salient composite index of biomarkers of stress and (5) to study the effect on stress of observance to the program at follow-up.

The following anthropometric measures will be obtained according to the recommendations of the International Society for the Advancement of Kinanthropometry[54]: standing height (m) and body mass (kg), waist circumference (cm), and lower limb bone length/breadth (cm).[54]

Body composition (muscle mass and fat mass) will be measured using DXA (DXA, QDR-4500A, Hologic, Inc., Waltham, MA, USA).6 and impedancemeter.80

Skin conductance will be measured using wristband electrodes with sampling rates at 2, 4, 8, 16 and 32 Hz during phases 1 to 3. The SC sensor (Q-Sensor®-Affectiva®, Massachusetts Institute of Technology, USA) is set on a wristband and has a 24-hour battery life when logging. In addition, it measures wrist movements with a built-in 3-axis accelerometer.

Blood flow velocity and myocardial longitudinal strain will be measured using speckle echocardiography (Vivid Q, GE Healthcare Biosciences, Piscataway, NJ, USA). All 2-dimensional (2D), time-motion, Doppler and 2D-strain acquisitions and measurements will be performed according to recent guidelines.[55,56] Left ventricular (LV) volumes and ejection fractions will be measured using Simpson’s biplane method.[57] LV mass will be calculated using the Devereux formula and indexed for height (Cornell adjustment). Pulsed Doppler LV transmitral velocities, including early and atrial waves, will be obtained using the apical 4-chamber view. Tissue Doppler imaging measures of myocardial systolic, early diastolic and atrial velocities will be assessed at the mitral annulus level in the apical 4- and 2-chamber views. The early transmitral flow velocity to the early diastolic tissue velocity ratio was used as an index of LV filling pressure. [58] Left atrium volume will be assessed using apical 4- and 2-chamber views. A graduation of LV diastolic dysfunction will be obtained according to recent guidelines. [58] 2D cine-loops (frame rate >70 ips) of at least five cycles will be recorded in the short-axis views (base, mid, apex), as well as in the apical 4-, 3- and 2-chamber views. 2D-strain analysis will be performed post-processing using EchoPAC 201TM software (GE Healthcare Biosciences, Piscataway, NJ, USA). Longitudinal and circumferential strains and strain rates, as well as apex and base rotations, will be directly obtained from the six-segment model. Twist mechanics will be computed from apical and basal rotational data using dedicated software (Scilab, Paris, France). For each view, the three cardiac cycles displaying the best image quality will be selected. Blood pressure and heart rate will be continuously monitored, and the systolic meridional wall stress, an index of afterload, will be calculated. LV end-diastolic volumes will also be obtained as a preload index.

Blood samples will be collected by a qualified paediatric nurse after participants have fasted overnight. Blood will be collected using a venepuncture at the brachial vein. After collection, blood will be centrifuged and aliquots will be stored (−80°) for subsequent analysis.

(triglycerides,6 cholesterol,6 low-density lipoprotein cholesterol,6 high-density lipoprotein cholesterol,6 glycaemia,6 insulin,6 ultra-sensitive C-reactive protein,6 cortisol,17,18 DHEAS17,18) as well as all other biochemical measures (leptin,6 brain-derived neurotrophic factor,20,87 interleukin-1b,6 interleukin-6,6 interleukin-1,6 tumour necrosis factor a,6 and neuropeptide Y19) will be assessed in the biochemistry laboratory of Clermont-Ferrand University Hospital. All analyses will be conducted by the same technician. Polymorphism of the angiotensin converting enzyme8-11 and polymorphism of serotonin12-14 will be measured from blood cells, as welll as telomere length.88

Questionnaire measures of the following variables will be obtained: stress and fatigue (visual analogue scale of 100 mm),[4] burnout (Maslach Burnout Inventory),[59] depression and anxiety (Hospital Anxiety and Depression Scale,[60] State–Trait Anxiety Inventory (STAI-Y),[61] 7-item Hamilton Anxiety Rating Scale),[62] mindfulness (Freiburg Mindfulness Inventory),[22] coping strategies (Brief COPE Questionnaire),[63] emotions (Emotion Regulation Questionnaire),[64] perception of work (Karasek’s Job Content Questionnaire),[65] self-efficacy (Perceived Self-efficacy Scale),[66] addiction to work (Work Addiction Risk test),[67] alexithymia (20-item Toronto Alexithymia Scale),[68] illness perception (Brief Illness Perception Questionnaire, adapted for stress at work),[69] metacognition (Metacognitions Questionnaire-30),[70] general health (12-item General Health Questionnaire (SF-12)),[71] lifestyle (e.g. smoking, alcohol consumption)[4] and physical activity (Recent Physical Activity Questionnaire).[72]

Bone microarchitecture will be measured using pQCT (XCT 3000 Stratec Medizintechnik Pforzheim, Germany).[73-75] Bone mineral content (BMC) (g/cm), volumetric cortical and trabecular BMC (mg/cm, total area (mm, cortical and trabecular area (mm and density (g/cm, and bone strength (mm will be assessed at the distal (4%) and proximal (66%) sites of the non-dominant tibia and radius. A planar scout scan was first conducted to determine the anatomical reference line for both the radius and tibia. Tomographic slices of 2 mm thickness were obtained at the 4% and 66% sites measured distally. Scan speed and voxel size were 20 mm/s and 0.4 mm, respectively. To ensure measurement quality, calibration checks will be performed by scanning a standard phantom with known densities prior to each scan. To calculate BMC, volumetric cortical BMC, volumetric trabecular BMC, cortical bone area, cortical density, trabecular area, trabecular density and stress strain index will be analysed using the Stratec pQCT manufacturer’s software. Contour mode 1 with a threshold of 180 mg/cm3 will be used to separate soft tissue and bone to analyse trabecular bone. Cortical bone will be identified and removed using a constant default threshold of 710 mg/cm3. A contour mode 3 with peel mode 1 at a threshold of 40 mg/cm3 will be used to assess muscle and fat cross-sectional areas.

Bone densitometry will be conducted using DXA (QDR-4500A, Hologic, Inc., Waltham, MA, USA). Bone mineral density (g/cm2), BMC (g) and bone area (cm2) will be determined for each participant. The DXA measurements will be taken for the whole body, the lumbar spine (L2–L4) and the non-dominant hip (including the femoral neck, and the trochanteric and intertrochanteric regions). All DXA scans will be conducted by the same technician and quality assurance checks will be performed routinely. The in-vivo coefficient of variation is 0.5%.

Statistical analysis will be performed using Stata software (version 13; StataCorp LP, College Station, TX, USA). All statistical tests will be two-sided and p < 0.05 will be considered significant. After testing for normal distribution (Shapiro–Wilk test), data will be treated either by parametric or non-parametric analyses according to statistical assumptions.

Mixed models will be used to analyse longitudinal data for fixed effects (before and after the residential thermal program), time-point evaluation and group × time interaction taking into account between- and within-participant variability (as random effects). A Sidak’s type I error correction will be applied to take into account multiple comparisons. The normality of residuals will be examined using the Shapiro–Wilk test. If necessary, a logarithmic transformation will be used to achieve normality for the dependent outcome. Multivariable analyses will be carried out with adjustment for covariates fixed according to epidemiological relevance and physical activity. Particular attention will be paid to the covariates ‘time between inclusion and beginning of the intervention’ (variable for each subject) and ‘time between beginning of the intervention and the follow-up’ (fixed for each subject).

Analysis of non-repeated data will be performed using analysis of variance (ANOVA) or Kruskal–Wallis tests. When appropriate (omnibus p-value less than 0.05), a post-hoc test for multiple comparisons will be used: the Tukey–Kramer test following ANOVA and the Dunn test following the Kruskal–Wallis test. Categorical variables will be compared using chi-squared or Fisher’s test. The Marascuilo procedure will be performed for multiple comparisons. Relations between quantitative outcomes will be analysed using correlation coefficients (Pearson or Spearman, according to statistical distribution). Fisher’s z transformation and William’s T2 statistic will be used to compare correlations between variables and within a single group of subjects.

A sensitivity analysis will be used to study the statistical nature of missing data (random or non-random); that is, baseline characteristics of participants with complete follow-up and those lost to follow-up will be compared using the aforementioned statistical tests. Even if random-effects models permit the analysis of data when a participant is lost to follow-up, the most appropriate imputation data method according to the statistical nature of missing data (multiple imputation data, last observation carried out (LOCF)) will be applied. More precisely, if a participant drops out of the study before it ends, then his or her last observed score on the dependent variable is used for all subsequent (i.e. missing) observation points if the LOCF imputation approach is used.

Data will be stored in the principal investigator’s office on a password-protected computer only accessible to members of our research team. In the electronic database, participants’ names will be replaced with numeric identity codes. Blood samples will also be labelled with the numeric identity codes and samples will be stored in our laboratories. Only aggregate results will be reported, ensuring participants’ anonymity.

Both DXA and pQCT provide measures of body composition and bone properties by exposing participants to low-level radiation: 0.0056 mSv from DXA scans (whole body, lumbar and hip) and 0.0014 mSv from pQCT scans (tibia and radius measures).[76] Over the duration of the study, the effective dose of 0.014 mSv will be administered.

The ThermStress protocol complies with the ethics guidelines for clinical research (Declaration of Helsinki) and has been approved by an ethics committee (CPP Tours Région centre Ouest 1, France - ANSM: 2016 A02082 49). The protocol has also been registered with the database ClinicalTrials.gov (NCT03536624). In accordance with ethical considerations, the chief investigator will be responsible for ensuring that participants understand the potential risks and benefits of study participation. Moreover, the chief investigator will be responsible for obtaining written consent from participants. The results will be disseminated at several research conferences and published in peer-reviewed journals.

Concluding remarks

The ThermStress protocol was designed to provide a better understanding of the effect of a short residential spa program to improve HRV and prevent work-related stress. The idea is to adapt a 3-week thermal cure for occupational burnout to a 6-day program, which is more accessible to workers and hopefully will produce similar results for work-related stress and burnout prevention. The long-term success of lifestyle interventions such as those proposed in work-related stress prevention depends on treatment observance (psychology, physical activity, nutrition).[77] We previously demonstrated that spa programs can play a major role in prompting sustainable lifestyle changes.[28] This shorter program may lead to better treatment observance during the 1-year follow-up, as psychological interventions, physical activity and diet will be supervised in the spa resort and participants will be accompanied by health care professionals. Findings from this protocol are expected to offer a more appropriate thermal program for burnout and work-related stress prevention that is accessible to the public.