**Assessment of the internal dynamics of thylakoids.** The result of neutron spin-echo spectroscopy measurements is the normalized intermediate scattering function *I* (*q*, *τ*) / *I* (*q*, 0) describing the pair-correlation distances observed in reciprocal space *q* = 2π/*d* and the time-dependent relaxation of such correlations due to the internal mobility of the sample investigated. As the first step in the analysis, the normalized intermediate scattering functions experimentally obtained, were fitted by a stretched exponential function with the stretching exponent *β* as a free fitting parameter. Figure 1 shows the *q* dependence of the stretching exponent. The values of *β* range between 0.2 and 0.4 in the small *q*-regime ≤ 0.07Å−1 with the rest of *β* around 0.6 for higher *q* values.

One can also observe the large *β*-errors, providing insights how stretched exponential function with a free exponent might not be the correct model to describe the experimental relaxation curves, especially for the data in the high *q*-regime. Therefore, in the next step of the evaluation we fixed the value of *β* = 2/3 according to Zilman & Granek18,19 model (abbreviated ZG from here on) proposed for bilayer membrane dynamics and the intermediate scattering function observed by NSE were fitted accordingly. ZG model provided a reasonable description of the NSE experimental data, with χ2 = 0.019 per degree of freedom, see Fig. 2.

The corresponding time information obtained is expressed as either relaxation time *τ* or inversely, as relaxation rates *Γ*ZG and calculated using Eq. 1, see Methods. In the following we used the spatial dependence of the relaxation rate *Γ*ZG to characterize the dynamics of our system (Fig. 3, both panels). The *Γ*ZG behavior as a function of *q* shows different patterns between the small *q*-range (black symbols) and the large *q*-range (red symbols) in Fig. 3, left panel. In the small *q*-range, there is no change in the relaxation rate values as a function of *q*. This points toward a relaxation process restricted by space, like diffusion in confined geometry. In the large *q*-range the relaxation rate is characterized a by a power law with the exponent of *m* = 3.88. A linear fit of all data yields a power law with the exponent of *m* = 3.28, when all *Γ*ZG values are used.

Given the value of the exponent close to *m* = 3, the *Γ*ZG /*q*3 behavior as a function of *q* is commonly represented to observe the linear dependence at high *q*-values, a typical signature of pure undulation motions in lipid bilayer membranes that helps extract the effective bending modulus of the membranes. It should come as no surprise that the dependence displayed in Fig. 3, right panel, shows a poor linearity between the relaxation rates *Γ*ZG /*q*3 at high *q*-values, ≥ 0.07Å−1. This points to additional relaxation due to superimposed dynamic processes on top of the simple undulation motions. The complex nature of the relaxation rate pattern arises from the complexity of the thylakoids stacked architecture within the chloroplast of the leaf and its intricate dynamical behavior. The photosynthetic membranes closely appressed together and enclosed in chloroplast feel the presence of other neighboring membranes in terms of fluctuation dynamics and one cannot easily deconvolute the single bilayer membrane undulation dynamics from the overall bilayer-bilayer interaction within the stack. The original ZG model is expected to be valid in the high *q*-regime, where 1/*q* < < interthylakoidal distance. The deviations observed in *Γ*ZG */q**3* vs *q* represents a limitation of the ZG model for characterizing bending fluctuation in complex architectural membrane structures, especially in living thylakoids where a multitude of other components beside bilayers are present and may contribute with additional relaxation modes.

**Characteristics of bending fluctuations.** The *q*3 dependence of the relaxation rate *Γ*ZG is controlled by the elastic modulus \(\stackrel{\sim}{\kappa }\), with higher bending coefficients corresponding to smaller relaxation rates18,20. The effective bending coefficient was calculated according to Eq. 2 from the intercept value *y*0 in Fig. 3, right panel (linear regression with slope = 0 depicted by the green dashed line) and expressed in units of thermal energy. Large value of \(\stackrel{\sim}{\kappa }\) >> *k*B*T* (\(\stackrel{\sim}{\kappa }\) = 1669.7 *k*B*T*) indicates a rigid system, a stack of membranes closely appressed together or membranes encapsulated in a tight structure with restricted space for free undulation motions, as well as it can indicate the presence of other structures like macromolecules grafted on the surface of the surfactant membranes that act like fluctuation suppressors18,20. All these are valid if we consider the architecture of thylakoids encapsuled by chloroplasts in the cells of *Lemna minor* leaf. The effective bending coefficient in stacked membranes depends also strongly on the viscosity of the solvent confined between membranes, as shown previously in numerous studies33. A rescaling of the bulk D2O viscosity, *η*D2O, to 3·*η*D2O as suggested by these studies brings the effective bending coefficient to \(\stackrel{\sim}{{\kappa }}\) = 185 *k*B*T* (Table 1), well within the extensive range observed for surfactant membranes, where values between 0.1–2000 *k*B*T* have been reported by various research.

Table 1

**Characteristics of bending fluctuations.** The bending modulus was calculated for *η*D2O = 0.00125 kg·m·s− 1 as D2O viscosity at 20°C, and for 3·*η*D2O as suggested by literature20,33.

Sample name | Γ/q3 (Å3/ns) | \(\stackrel{\sim}{\varvec{\kappa }}\) (kBT) | \(\stackrel{\sim}{\varvec{\kappa }}\) for 3 · *η*D2O (kBT) |

Lemna Minor | 1.98 ± 0.33 | 1669.7 ± 556.6 | 185.5 ± 61.8 |

Given the above, an analogy can be made between the observed excess dynamics of thylakoids and the excess dynamics related to shape fluctuations in oriented lamellar phase microemulsions34,35. Dynamics of oriented lamellar phases model implemented in the software Jscatte*r*36,37 calculates the coherent intermediate scattering function arising from dynamics of oriented lamellar phases at the length scale of the intermembrane distance (see Eq. 16 in reference Mihailescu *et al*., 200234). In our case the model provided information on the single membrane apparent bending modulus in thermal energy units, the membrane thickness, and solvent viscosity, using the film-to-film distance (241.65Å, first SANS peak @ *q* ~ 0.026Å−1) representing the periodicity of the structure as input parameter (Fig. 4).

The best fit was obtained assuming *q*-dependence of the bending modulus and solvent viscosity, with χ2 = 0.017 per degree of freedom (117 dof, with 31 parameters fitted simultaneously for the ten experimental scattering functions, *i.e. q*’s). The bending modulus calculated for single membrane varies between 0.59 *k*B*T* and 3.72 *k*B*T* (Fig. 5), higher but not too different from values determined for bi-continuous surfactant phases38 (~ 1.3–1.5 *k*B*T*). Single bilayer thickness calculated is 65 Å ± 24 Å and strongly corelated with the solvent viscosity. In the analysis only bulk apparent viscosity was assumed and fitted for each *q* to account for local viscosity fluctuations. These variations can be seen in Fig. 5. A mean value of 0.75 mPa*s (with a large standard deviation of 0.75) can be calculated from the viscosity values, falling well within the range observed for the n-alkanes and chloroalkanes, as well as aqueous solution with high solutes concentration39. We theorized that the sensitivity of the calculated model parameters to observable bulk viscosity represents a qualitative description of local viscosity fluctuations due to anchored proteins on the membrane surfaces, increased local friction with the solvent and variabilities in the membrane internal viscosity due to variations in lipids concentration.

**Characteristics of local fluctuations.** The spatial dependence of the decay rate *Γ* shows excess mobility in addition to the underlying *q*3-dependent undulation dynamics. To quantify these deviations we used the approach established by Nagao and collaborators13,21 where the excess motion is described as local shape and thickness fluctuations with peristaltic and protrusion motions, see Eq. 3 in Methods. The intent is to describe the trend of the relaxation rate by Lorentz functions, using the previous peak model analysis used for photosynthetic cells1 based on simulation results for thickness fluctuations21 and protein-induced deformations in lipid bilayers20,40. There are two strong deviations identifiable, one in the small to intermediate *q*-range where we sample distances proportional to the center-to-center bilayer distance, and another deviation in the high *q*-regime at distances in the range of a single membrane thickness. In the past, similar excess dynamics have been attributed to swelling of the membrane pair, *e.g.*, thickness fluctuations, bilayer-bilayer interactions, and dynamics due to protrusion and diffusion of proteins at the membrane surface. The decay rate *Γ* experimentally obtained fit was separated again, for the purpose of this study, in two *q*-regimes: a high *q*-regime ≥ 0.07Å−1 (red data in Fig. 6) and a small *q*-regime ≤ 0.07Å−1 (black data in Fig. 6) and two Lorentzian were used for the analysis of the two *q*-regimes. The best fit was obtained by Lorentzian functions having the peak center situated around the same position as the first and the last observable SANS correlation peaks: *q* ~ 0.027Å−1 and *q* ~ 0.11Å−1. This is a clear indication of excess dynamics occurring at length scales corresponding to both: bilayer periodicity and single membrane thickness. The calculated parameters from the Lorentz fit are presented in Table 2.

The large difference observed in the decay rate of the local fluctuations, *A*·*q*03, between small and high *q*-regime is a good indication of the underlaying dynamical process. In the small *q*-regime the faster decay rate points toward longer relaxation times and slower motions. These slow motions can be attributed to progressive swelling of the thylakoid pair stack as well and to adjacent bilayer repulsive interactions. The photosynthetic membrane pair is a living structure, and a certain amount of flexibility is needed to adjust to the active photosynthetic process. The swelling of the thylakoid pair happens within the restricted space of the chloroplast, up to the point where is contracted by the adjacent bilayer swelling. Although we lack a more discrete refinement of *Γ vs q* dependence, as well as values at the very low *q*, we can still formulate some hypothesis on a reasonable behavior for living photosynthetic membranes under restricted spatial architecture within chloroplasts.

Table 2

**Characteristics of local thickness fluctuations.** *A* = damping frequency of the peristaltic mode, *ξ* = the peristaltic mode amplitude, *q*0 = the peak maximum position, local length scale of thickness fluctuations. The product *A*·*q*03 describes the decay rate of the local fluctuations. **Note that the parameters are calculated for each peak displayed in* Fig. 6 *using the intercept value y**0* *= Γ*ZG*/q*3 *= 0.89 of the two Lorentzian.*

Sample name | ΓZG/q3 (Å3/ns) | q0 (Å−1) | A (Å3/ns) | ξ (Å) | A*q03 x 10− 4 (ns− 1) | χ2 |

Peak1 | 0.89 ± 0.02 | 0.027 | 0.24 ± 0.002 | 0.03 ± 2.5E-4 | 0.045 ± 0.0004 | 3.4E-4 |

Peak2 | 0.89 ± 0.02 | 0.109 | 0.13 ± 0.03 | 0.04 ± 0.01 | 1.74 ± 0.035 | 0.42 |

With increasing distance from the marked high *q* = 0.11Å−1 (*d* ≈ 50Å) to the marked low *q* = 0.027Å−1 (*d* ≈ 250Å), the field of view changes from the local single membrane leaflet to center-to-center distance of adjacent bilayer membrane-pair stack. The thermal fluctuations can induce collisions between neighboring membranes. These collisions give rise to repulsive bilayer-bilayer interactions caused by the reduction of entropy, called Helfrich steric repulsion41. Experimental work on lamellar microemulsion of sodium dodecyl sulfate (SDS) of Safinya *et al*.,42 showed that for intermembrane distances between ~ 40Å and ~ 170Å Helfrich steric interaction is the dominant interaction. The *q-*window that NSE was able to sample during the measurements on *Lemna minor* 0.04Å−1 < *q* < 0.13Å−1 with corresponding distances of 157Å < *d* < 48Å sits within the exact region of steric repulsion described above. In the high *q*-regime, however, the larger decay rate *A*·*q*03 (38-fold faster), points toward much faster motions. The sole existence of an observable deviation at higher q (distances in the range of a single membrane thickness, ~ 50Å), points toward short-range dynamics due perhaps to protrusion and diffusion of large proteins anchored on the surface of the single leaflet membrane, that contaminate the observed NSE relaxation time window. Future planned experiments at extended *q*-range with the characteristic dynamical trends combined with a thorough analysis of the SANS correlations peaks under variable illumination conditions are planned and will expand further our understating of these photosynthetic membranes dynamics over an extended *q*-scale.