This single-center prospective study was approved by our institution’s ethics committee (CE 19.358) and was registered to Clinicaltrials.gov (NCT04716101. Registered 14 January 2021. https://clinicaltrials.gov/study/NCT04716101). Written informed consent was obtained from all participants.

Between February 2021 and June 2022, individuals meeting inclusion criteria—adults aged 18–75 with low back or referred pain above or just below the gluteal fold (≥ 3/10 on a numerical pain rating scale) lasting at least 3 months, present at least 50% of the time during the day—were recruited by fellowship-trained anesthesiologist (A.B.) or musculoskeletal radiologist (N.J.B.) with 25 years of experience, at our institution’s Pain Clinic and Radiology Spinal Intervention Unit. Exclusion criteria included previous back or lower extremity injury, spinal surgery, pain attributable to specific known pathology, corticosteroid use or lumbar spine injection in the past 3 months, and pregnancy. Volunteers without a history of low back pain or chronic pain limiting daily activities and meeting exclusion criteria were recruited through advertisements.

A research assistant collected demographic data and measured weight, height, and lumbar spine range of motion following the Schober test method by Tousignant et al [16]. Additionally, participants completed three self-administered questionnaires: the International Physical Activity Questionnaire [17], the Brief Pain Inventory [18], and the Oswestry disability index [19].

The radiologist conducted radiofrequency (RF) data acquisitions using a system equipped with a 12L5 linear array transducer (Terason t3000 v4.7.7, Terason Ultrasound, Burlington, USA). Participants were prone on an articulated motorized table (Echo Flex Model 4800, Ibiom Instruments, Sherbrooke, Canada), aligning the iliac crests with the hinge point of the table (Fig. 2).

Ultrasound scanning. a A urethane resin plate, one cm-thick and rectangular (depicted by the black rectangle), with a customized hole to accommodate the transducer’s surface, was positioned on the paraspinal muscles at the peak point of the ES muscle, specifically at the L2–L3 interspinous level. Echo gel (blue layer) was used, and the radiologist conducted sagittal plane scans while utilizing the plate to stabilize the transducer. b Participants lay prone on an articulated motorized table, aligning their iliac crests with the hinge point. Scans were performed while the lower extremities underwent passive movement, as the distal portion of the table descended 20° before returning upward to the neutral position

At the L2–L3 interspinous level, we identified the peak point of the ES muscle on a transverse B-mode image and marked it on the skin (Fig. 3). Subsequently, we placed a one cm-thick rectangular-shaped urethane resin plate featuring a hole tailored to the transducer’s surface at the marked position. The hole was filled with echo gel. The radiologist used the plate as a stabilizer to prevent longitudinal and lateral transducer movements during scanning. The image (range 3.0–6.0 cm) and focus (range 1.3–4.5 cm) depths were adjusted to the participants’ body habitus with frequency preset “middle” and frame rate between 30.6 Hz and 48.6 Hz.

Ultrasound imaging and layer identification of paraspinal soft tissues. a Transverse (short-axis) and (b) longitudinal (long-axis) B-mode images of the right paraspinal soft tissues in a 50-year-old man with NSLBP. a The transverse image at the L2–L3 interspinous space (large arrow) shows the ES muscles with their hyperechoic aponeurosis (thin arrow) and the overlying TLF (open arrow). The peak convexity of the ES muscles (asterisk) was marked on the skin, and scanning was performed at this marked position in the longitudinal plane. b The longitudinal image displays the paraspinal tissues’ ultrasound anatomy, including the dermis, subcutaneous adipose tissue, multilaminar TLF, hypoechoic hyaluronan-rich loose connective tissue layer, hyperechoic ES aponeurosis, and muscles. For outlining paraspinal tissue layers, the medical students placed dots along the soft tissue layer margins on the image center’s first frame of the table-stationary phase. An algorithm then measured the thickness of the subcutaneous adipose tissue, TLF, and juxtamuscular zone (vertical green lines). The juxtamuscular zone included the hypoechoic hyaluronan-rich loose connective tissue and the hyperechoic ES aponeurosis. Distances were recorded in pixels and converted to millimeters using a 0.026 mm/pixel resolution. On the middle frame of the table-downward phase, two 10-mm red lines were positioned at the upper TLF and lower ES aponeurosis margins (shown together for convenience), defining the region of interest for ShS calculation. TLF, thoracolumbar fascia; NSLBP, nonspecific low back pain

The radiologist scanned in the longitudinal plane, adjusting the transducer perpendicular to the hyperechoic ES aponeurosis. During passive movement of the lower extremities, a 10-s cine loop was captured as the distal part of the table descended 20° at an angular speed of 7°/s, followed by an upward return to the neutral position. This 10-s cine loop comprised a 4-s table-stationary phase, a 3-s table-downward phase, and a 3-s table-upward phase. Both sides of the spine were scanned three times, totaling six cine loops per participant, with each acquisition spaced one minute apart.

RF datasets were analyzed with custom MATLAB programs (vR2022a, Mathworks Inc., MA, USA). The radiologist trained two medical students (S.D., A.A.) to outline anatomical structures on B-mode reconstructed images from ten consecutive US cine loops. All were blinded to the participant group. Subsequently, the students completed the same task for the remaining datasets. On the first frame of the table-stationary phase, they measured the thickness of the subcutaneous adipose tissue, TLF, and juxtamuscular zone encompassing the loose connective tissue layer and the ES aponeurosis. Additionally, on the middle frame of the table-downward phase, they traced a 10-mm long horizontal line at the superficial margin of the TLF and the deep margin of the ES aponeurosis. The region between these two lines defined the region of interest (ROI) used for ShS calculation (Fig. 3).

The TLF’s ShS was computed using the Lagrangian Speckle Model Estimator, an ultrasound elastography analysis algorithm. This method was previously validated [20] and used to assess atheromatous plaque vulnerability [21,22,23], lung kinetics [24], and single-cell mechanics [25]. Unlike methods relying on B-mode images, this algorithm calculates parameters based on RF data, ensuring that ShS results are unaffected by scanner settings and image post-processing techniques, thus improving accuracy. Additionally, RF processing enhances the precision of tissue movement tracking.

The analysis involved the measurement window sliding over the ROI in the table-downward and upward phases of RF data. A cross-correlation method compensated for translating the measurement window between successive RF images. Subsequently, 2D strain components were computed using an extended version of the optical flow equation [20]. These components included axial and lateral translations, strain, and ShS, where “axial” referred to the direction along the ultrasound beam (vertical in the image), and “lateral” denoted the horizontal direction perpendicular to the axial direction in the image and the anatomical longitudinal direction [23, 26]. We extracted the time-varying instantaneous absolute lateral ShS among those parameters to characterize the TLF’s sliding movement. Here, “instantaneous” denotes measurements between successive RF frames, and the absolute value signifies the magnitude of the shear movement, irrespective of its direction.

Vibrations at the start and end of the table’s movements caused axial translation noise. We used a time-varying axial translation curve to identify stable intervals of lateral TLF movement, excluding peak translations (Fig. 4). The resulting intervals in the table-downward phase [NSLBP 100.49 ± 13.87 vs controls 105.78 ± 16.26 (p = 0.27)] and the table-upward phase [NSLBP 98.45 ± 16.51 vs controls 99.94 ± 17.32 (p = 0.41)] had a similar number of frames. We computed the cumulative absolute lateral ShS magnitude (C|ShS|L) and the maximum absolute lateral ShS (Max|ShS|L) over these intervals.

Axial translation curve over time during table downward and upward movements. This graph shows the absolute vertical displacement over time for a 74-year-old asymptomatic male volunteer, highlighting the stationary, downward, and upward table phases. Note that only a segment of the 10-s cine loop is shown. The arrows indicate that peaks in axial translation are due to table vibrations at the start of the downward movement, at the table reversal point, and at the end of the upward movement. For RF data analysis, the stable intervals during the downward and upward table movements were determined by excluding these three axial displacement peaks. Consequently, the table’s downward and upward intervals are not continuous and are analyzed separately

Intra-operator reliability for data acquisition and computation was evaluated using table-downward interval RF data from the pre-and post-sham intervention cases, encompassing data from 15 NSLBP and 16 control participants.

Participants were randomly assigned in a 1:1 ratio to receive either a standardized massage therapy or a sham technique. A certified therapist (N.T.) with 10 years of experience performed the massage technique by applying deep and slow compression and shearing forces to the participant’s lower back, moving from the iliac crest to the lower ribs, following inter-jurisdictional practice competencies guidelines [27]. This included ten strokes along the right and left ES muscles’ medial, middle, and lateral parts, totaling 30 strokes per side. A custom-made force sensing system comprising a unidirectional reusable flexible force sensor (FlexiForce Standard Model A301, Tekscan, MA, USA) affixed to a massage tool (Deep Pressure Thumb Saver Massager, MIGLEO, Shenzhen, China) was used to measure and calibrate the force during the massage strokes. Connected to a computer through an amplifier (non-inverting op-amp circuit) and an analog-to-digital converter (NIUSB-6251, National Instruments Corp. TX, USA), the force sensor gauged the therapist’s thumb’s force on the massage tool in Newtons using a pre-calibrated force-voltage relationship. Before the massage, the therapist determined the pressure threshold causing discomfort. During the massage, the therapist exerted pressure between the upper discomfort threshold and two-thirds of that value to ensure participant comfort. We standardized the force applied to relative intensity, mirroring clinical practice. The therapist performed the sham technique by placing one hand on the ES muscles without exerting pressure sequentially on the right and left ES muscles. Both interventions were conducted over a hospital gown, each lasting 2.5 min. Ultrasound data acquisitions were repeated immediately after the intervention.

The study included a convenient sample of 60 participants equally divided into two groups, with 32 participants per group to account for potential failed examinations. Descriptive statistics were used to analyze demographic and clinical characteristics, ShS parameters, and the thickness of the paraspinal soft tissues. The Shapiro–Wilk test assessed data distribution, and square root transformations were utilized for non-normally distributed data. Linear mixed-effects models were employed to analyze between-group comparisons of ShS parameters, examine the impact of the interventions on ShS parameters, and investigate the correlation between ShS parameters and pain and disability scores. Intra-operator reliability was assessed using the intraclass coefficient correlation (ICC) with 95% confidence intervals. One author (MHRC) conducted the statistical analyses using R software (v4.2.1, The R Foundation for Statistical Computing, Vienna, Austria) with a two-tailed test and a significance level of 0.05.