Imaging and characterization of transitions in biofilm morphology via anomalous diffusion following environmental perturbation

: Microorganisms form macroscopic structures for the purpose of environmental adaptation. Sudden environmental perturbations induce dynamics that cause bacterial biofilm morphology to transit to another equilibrium state, thought to be related to anomalous diffusion processes. Here, detecting the super-diffusion characteristics would offer a long-sought goal for a rapid detection method of biofilm phenotypes based on their dynamics, such as growth or dispersal. In this paper, phase-sensitive Doppler optical coherence tomography (OCT) and dynamic light scattering (DLS) are combined to demonstrate wide field-of-view and label-free internal dynamic imaging of biofilms. The probability density functions (PDFs) of phase displacement of the backscattered light and the dynamic characteristics of the PDFs are estimated by a simplified mixed Cauchy and Gaussian model. This model can quantify the super-diffusion state and estimate the dynamic characteristics and macroscopic responses in biofilms that may further describe dispersion and growth in biofilm models.


Introduction
Distinct features of bacteria in biofilms, compared to a single bacterium, arise via denser cellular communications [1][2][3]. The collection of cells increases the concentration of communication molecules and the collective cellular responses become more systematic in contrast to single-cell responses [4]. Most bacterial strains favor forming a biofilm to develop environmental resistance by embedding themselves within an extracellular matrix and minimizing metabolism [5,6]. The biofilm structure mitigates the susceptibilities of bacteria to antibiotics and immune responses that lead to a chronic infection in the host [7,8]. Such microbial resistance can complicate treatment and recovery from acute and chronic infections, and in rare cases cause lethal health outcomes to immunocompromised and elderly patients. Various bactericidal and biofilm formation inhibition methods have been proposed, but the frequent genetic mutations and dormant bacteria in biofilms often mitigate their effectiveness [9][10][11][12].
To understand and control the mechanisms of biofilm formation mediated by cellular communications, quorum sensing (QS)-driven molecular signaling for biofilm formation and disintegration has been investigated and previously demonstrated [13][14][15]. Biofilms have a mechanically stable structure and bacterial growth or dispersal dynamics induce structural deformation [16]. Quantifying these dynamic properties would provide deeper understanding of bacterial behaviors and also potential applications for controlling biofilms. However, the conventional observations of biofilm dynamics by measuring biofilm thickness and distribution often take hours or days [17,18] and the development of a more instantaneous measurement method would be preferred. Biofilm formation and dispersion also display unique dynamic features in response to external perturbations [19][20][21]. Various dynamics occurring in biofilms after chemical stimulation [22][23][24][25][26] can be investigated by a random walk model interacting with complex geometries, which are slowly varying by environmental perturbations [27][28][29], and a long relaxation time of the expanding or dispersal dynamics may result in anomalous diffusion [30,31]. Not only measuring dynamic responses of biofilm, but also seeking possible correlations between quantitative properties of biofilms and fundamental dynamic processes has been a holy grail in the biophysics community [32], and such characterization may be feasible by non-invasive and label-free optical measurements [33][34][35][36].
By combining dynamic light scattering (DLS) and phase-sensitive Doppler optical coherence tomography (OCT) techniques, it is possible to acquire high-resolution dynamic images of biofilms without the use of any perturbative labels or dyes [37][38][39][40][41][42][43][44][45]. The mean and variance of temporal phase displacements acquired from the dynamic images can be correlated to the mean and variance of particle displacements [46]. When perturbed, the probability density functions (PDFs) of biological dynamics measured by Doppler shift show anomalous diffusion characteristics, and the alternations of the anomalous diffusion characteristics can provide important information to understand the internal dynamics of the sample [35,36]. Many biological processes have anomalous random walks whose displacement PDFs have power-law tails caused by rare ballistic flights [47], and such anomalous random walk characteristics can be quantitatively estimated by acquiring the logarithmic slope α of the power-law tail [48]. Anomalous diffusion can be induced by many physical factors such as an external force field, active transport, or interactions with complex geometries that result in the statistics that deviate from the Markovian process and generate Lévy walk [31]. The physical interpretation of such behaviors is called a memory effect that causes a longer correlation of stochastic processes [30,31]. For instance, a random walker of a Gaussian process always has independent stochastic paths, however, in contrast, a random walker under an external driving force field generates biased and correlated paths. Measurements of the anomalous characteristics could provide insight into the fundamental properties of system dynamics but also demonstrate dynamic phenotypes for biomedical applications [35,36].
Biofilm growth and the associated dynamics is a process of adaptation to the specific environment [18,49,50]. This process results in different biofilm morphologies when grown on an air-agar interface or in a liquid medium [51]. Alterations in the environment introduce transitions in biofilm morphologies that result in dispersal and expansion. Such transition processes may have anomalous random walks [29,52,53] by a sufficiently large relaxation time and the memory effect [31]. To be more specific, when a biofilm grown on an air-agar interface is submerged in a medium, the biofilm slowly swells and the tightly embedded biofilm structure consisting of a polysaccharide extracellular matrix will slowly disperse. The biofilm will release random-walking dynamic particles that used to be a part of biofilm during the process. The expected schematic diagram of a biofilm grown on an air-agar interface and submerged in a medium is shown in Fig. 1(a). The structural and morphological changes of the biofilm after being submerged in a medium are shown in Fig. 1(b), and more detailed dynamic characteristics were acquired by the phase-sensitive OCT system in Fig. 1(c).
The volumetric deformations of the submerged biofilm ( Fig. 1(d)) and mediating dynamics of the random walkers exerting structural deformation were investigated by detecting the Doppler shift that occurred by dynamic light scattering ( Fig. 1(e)) [35,36]. The biofilm in a medium is assumed to release dynamic particles behaving like random walkers that can freely move and be trapped again while interacting within the biofilm. The expected trajectory of a dynamic particle has different characteristics when it is traveling in free space versus when confined or trapped within a biofilm. Such trajectories can be modeled by Lévy walks that have localized diffusion and rare ballistic walks at the same time [47]. The distributions of the phase displacements ∆ϕ of Characterizing dynamics within biofilms. (a) Schematic diagram of biofilm and dynamic particle trajectory (black solid lines) when exposed to a medium. Dynamic particles (yellow spheres) interacting with biofilm (green) can be either trapped or escape from the biofilm. While dynamic particles are interacting within a biofilm, the random-walk characteristics may change. The trajectory is a Lévy walk generated when α = 1.5 in Eq. (3). (b) Scanning electron microscope (SEM) and false colored 1X optical en face camera images of a S. aureus biofilm. In the SEM image, embedded bacteria formed a dense structure. The optical images show the biofilm (yellow and green) and agar surface (blue) before (0 min) and after (60 min) being submerged in a medium. The red dashed lines in the optical images represent the locations of OCT B-scans (2D scan) region. (c) Phase-sensitive Doppler OCT configuration, and (d) volumetric images of the biofilm before and after exposure to a medium. The applied medium induced structural deformation of the biofilm structure. Scale bars in (d) represent 1 mm, 1 mm, and 400 µm along the x, y, and z-axes, respectively. The relative height of the biofilm from the agar surface is represented in a false-color scale (blue at the agar surface). (e) Probability density functions (PDFs) of phase displacements of backscattered light were experimentally acquired from a biofilm submerged in a medium, and fitted by Lévy distribution. The PDF of the background noise determined by the intensity fluctuations shows a Gaussian-like distribution (α = 1.92). The PDFs of biofilm exposed to a medium show the power-law distributions where the dynamic biofilm is close to Cauchy distribution (α = 1). The estimation of the power-law tail slope of ∆ϕ −α−1 is drawn as a dashed line above the noise floor (flat solid line).
Doppler shift are shown in Fig. 1(e), where the phase displacements caused by the background noise are close to a Gaussian distribution and the dynamic biofilm is close to a Cauchy distribution. The dynamics of submerged biofilms, in general, show a super-diffusive Lévy distribution which is a distribution between Cauchy and Gaussian distributions. The PDF shapes of the power-law tails in the log-log plot show that each distribution exhibits orders-of-magnitude differences of the likelihood at a larger displacement as the power-law tail decays as ∆ϕ −α−1 . The accuracy of the power-law tail shape estimation is determined by the data size. As the log-log plot requires an exponentially increasing data size for accurate estimations of the shape of power-law tails, high-numerical-aperture optical configurations with a finite field-of-view may not provide a sufficient amount of dynamic information within a reasonable amount of sampling time [54]. The phase-sensitive OCT technique, however, is applicable, as data is acquired over a wide area with a fast sampling rate, and with a relatively high spatial resolution.
In this study, the dynamics of Staphylococcus aureus (S. aureus, non-motile gram-positive bacteria) biofilm induced by various reagents were measured by phase-sensitive Doppler OCT.
The bacterial distribution in biofilm shows a packed distribution embedded within an extracellular matrix shown in Fig. 1(b) [55], and the sudden exposure to a medium is assumed to generate dynamic particles. To exclude any inherent bacterial motility and to isolate the diffusion effect of dynamic particles, a non-motile bacterial species was selected. The conventional dynamic parameters of phase-sensitive Doppler OCT were acquired to demonstrate the relationships between the dynamic parameters and anomalous diffusion, experimentally, and to characterize the dynamic features of these biofilms. The PDFs of phase displacements induced by the Doppler shift were acquired and were fitted to the Lévy distribution with a fractional characteristic exponent function to estimate the stable distribution shapes between Cauchy and Gaussian distributions [48]. Also, the interpretations of the alternations of the power-law tails of the PDFs will be discussed.

Dynamic light scattering and phase-sensitive Doppler OCT
Light scattered from dynamic particles introduces a Doppler shift (ω D ) which carries the dynamic information [33,36]. When the light of a central wavelength λ 0 is illuminating a target along the z-axis within a refracting material of refractive index n, the relationship between the velocity component of the particle along the optical axis (v z ) and the induced Doppler shift ω D is where q is a photon momentum transfer vector and θ is a scattering angle. From this relation, the induced Doppler shift of backscattered light (θ =π, λ 0 =1325 nm) from a dynamic particle in water (n = 1.33) with a velocity of 1 µm/s along the optical axis (z-axis) is about 2 Hz (∼12.6 rad/s). Such minute frequency modulation induces phase modulation that results in temporally fluctuating beating interference. Practically, the phase value of a single voxel is determined by the vector sum of backscattered photons from multiple dynamic particles in a voxel, and the phase displacement between consecutive frames represents the average Doppler shift induced by the random walking dynamic particles [56,57]. The phase variance acquired from the time series of the Doppler shift is a key parameter of Doppler OCT that represents how motile the object is [33,36]. A higher variance value represents a diffusion process with a larger diffusion coefficient that is induced by more fluctuation of the phase angle of the backscattered photon [56,57]. The velocity estimated from the phase displacement of a voxel by Eq. (1) is defined as the Doppler velocity. The diffusion characteristics of dynamic particles are estimated by analyzing the shape of the phase-displacement PDF of consecutive frames [35,36]. The Doppler shift of the backscattered light is assumed to be the result of a single backscattering. The light-matter interaction with the scattering angle θ ∼ 0 would not induce a significant amount of Doppler shift by Eq. (1). The multiple back-and-forth photon scattering events were not considered because multiple mirror images were not observed during measurements [58]. The phase-sensitive Doppler OCT system in Fig. 1(c) is based on a spectral-domain OCT configuration which illuminates a target with low-coherence light from a superluminescent diode (λ = 1325 ± 50 nm, SLD135 Thorlabs), collects backscattered light by an achromatic lens (Thorlabs, AC300-050-C), and detects the low-coherence spectral interference by a spectrometer (SU-LDH2 1024-pixel InGaAs line-scan camera, Goodrich) [34]. The lateral and axial resolutions of the OCT system in free space are 8 µm and 16 µm, respectively. To acquire a stable signalto-noise ratio, data were oversampled and the real dimension of a pixel is approximately 4 µm by 4 µm along the lateral and longitudinal axis whose voxel size corresponds to 128 µm 3 . The complex-valued fields (E(t)) of OCT images were reconstructed from the interference in the spectral domain by the Fourier transform. To avoid phase wrapping while acquiring phase displacements (∆ϕ), consecutive frames (sampling period τ) were divided, and used to determine the phase angle closest to zero, as in The phase displacement represents the mean phase shift, represented by the Doppler shift over the sampling period τ. Note that the phase displacement corresponds to an average displacement of particles in a voxel, as the OCT signal is processed by collecting backscattered photons from particles in a voxel, and averaging their phase-angle contributions. From Eq. (1) and (2), the time series of ∆ϕ (radian) were acquired, and the mean and standard deviation of phase displacements were calculated. These parameters characterize the dynamic properties of the Gaussian diffusion. Practically, biological systems have various random processes with diverse dynamic parameters such as persistent length and velocities; therefore a single Gaussian estimation may be insufficient to fully demonstrate such information [35,36,59]. The characteristic parameters of the PDF were acquired by fitting the Lévy distribution [48] The Lévy α values of super-diffusion are between 1 and 2, which are between Cauchy (α = 1) and Gaussian (α = 2) distributions. These describe the distribution of displacements by Lévy walk whose trajectories have localized diffusion and occasional ballistic jumps. The parameter γ is a scaling parameter that determines the width of the distribution. The phase value of a coherently reconstructed image is the average phase angle of backscattered photons from the dynamic particles in a voxel within an exposure time (11 µs). The net Doppler shift induced by dynamic particles is measured by subtracting phase angles between frames. The subtraction rejects the constant background phase and isolates dynamic information. The Doppler velocity estimated by the phase displacement is given by The frame rate in this study was 55 Hz (1000 A-scans per frame) and the phase domain is -π to π. The measurable maximum amplitude of the Doppler velocity with a 55 Hz frame rate is about 13.7 µm/s. A higher speed that induces a phase shift beyond the phase domain would wash out interference and appear at the phase noise floor.

Phase jitter noise estimation
Optical sources have intensity fluctuations. Configuring an interferometer with the fluctuation source, the estimation of the phase displacement ∆ϕ in Eq. (2) has uncertainty as we are measuring the real part of the interference fringes. The acquisition of a PDF of ∆ϕ is sensitive to the phase noise floor [60][61][62]. The temporal intensity fluctuation creates the jitter noise of the phase displacements as the phase information is acquired from the spectral interference between the reference and object/sample arms. The estimated phase uncertainty between the consecutive frames is given by [63] whereĪ is the temporal and spectral mean intensity, and δI(τ) and δϕ(τ) are the temporal intensity and phase fluctuations, respectively, within the sampling period τ. To measure the phase-noise floor of the OCT system, a static phantom (silicone) target was scanned and the corresponding PDFs with different sampling periods τ were acquired. Frame rates were varied from 1.5 Hz to 76 Hz and the phase noise PDFs showed Gaussian-like distributions as the intensity noise had a Gaussian distribution. The PDF by the background noise showed α = 1.92 in Fig. 1(e) with the frame rate at 55 Hz, which is close to a Gaussian distribution. All data acquisitions were performed at this minimal-noise sampling frequency. As the statistics of the background noise have a Gaussian distribution, precise estimation of the subdiffusive diffusion (α ≤ 1) is difficult, so only super-diffusive diffusion (1 ≤ α ≤ 2) characteristics were considered in this study.

Lévy distribution estimation
The numerical estimation of the Lévy distribution was conducted by the maximum likelihood estimator (MLE) of the log-likelihood of the model function given in Eq. (3) for a careful estimation of power-law distributions [64]. The log-likelihood function l α,γ of the experimentally acquired PDF (p(∆ϕ)) and the model function L α,γ is defined as and was used to numerically acquire the optimized α and γ parameters of the model function L α,γ when the log-likelihood function l α,γ had a maximum value. Note that the experimentally acquired PDFs were not always continuous functions. Therefore, the numerical estimation was done in a discrete manner and Eq. (6) is denoted in a discrete form. Fitting results are shown in Fig. 1(e). The experimentally acquired PDF of the phase displacement distribution has symmetric and power-law tailed distributions (∼∆ϕ −α−1 ). To increase the accuracy of the MLE, the statistics of the phase displacements were collected from an adequate number of voxels and with a sufficiently long measurement time. The number of data points acquired from the scanning of the cross-sectional area of a biofilm is about 10 7 from 10 5 voxels (the average size of the cross-sectional area of biofilm ∼ 1000 by 100 voxels, corresponding to 4 mm by 400 µm along the lateral and depth axis) and 100 repeated temporal measurements at a 55 Hz frame rate. From the background noise measurement in Fig. 1(e), the effective dynamic range of PDF above the phase noise background is approximately 10 3 .

Biofilm preparation and perturbation induced dynamics measurement
Planktonic Staphylococcus aureus (S. aureus, NR51163, BEI RESOURCES) was cultured in a brain heart infusion medium (BHI, Difco TM , 7.7%, beef heart infusion 9.8%, proteose peptone 10%, dextrose 2%, sodium chloride 5%, and disodium phosphate 2.5%) for 24 hours and diluted in phosphate buffered saline (PBS, BP2944-100, Fisher BioReagents TM , 0.01 M phosphate buffer, 0.0027 M potassium chloride, and 0.137 M sodium chloride) by 10 −6 . A 100 µL volume of the diluted culture was widely spread on a BHI agar plate (Difco TM , 7.7%, beef heart infusion 9.8%, proteose peptone 10%, dextrose 2%, sodium chloride 5%, disodium phosphate 2.5%, and 15% agar) and incubated for 24 hours at 37 degree Celsius with 5% carbon dioxide (CO 2 ) and 70% humidity. When colonies formed on an agar plate, a whole colony was carefully collected and transferred to a new BHI agar plate. The agar plate was incubated for 23 ± 6 days at the same temperature, CO 2 concentration, and humidity. The maturation of biofilms was determined by the growth of the biofilm radius on the agar plate (thickness ∼ 300 µm, diameter ∼ 2 cm). The accessory gene regulator QS of S. aureus is known to trigger biofilm growth or the disperse phase that depends on the glucose level of the environment [17,65]. Glucose-rich and depleted environments introduce systematic growth and dispersal of S. aureus biofilm, respectively. Four different media were applied to S. aureus biofilms prepared on agar plates. To investigate the diffusion patterns of this bacteria in biofilms, phosphate-buffered saline (PBS), and brain-heart infusion (BHI) media, with and without the antibiotic (Abx, amoxicillin 10 µg/ml, A8523, Sigma-Aldrich) were used in this study [17,65]. A 5 ml volume of medium was added to the biofilm growing on BHI agar and 100 2D OCT B-scans (2-dimensional cross-sectional scanlateral vs. depth) were acquired at the same cross-sectional location at a scanning rate of 55 Hz. Each B-scan image consisted of 1000 A-scans (1-dimensional depth scan) with an exposure time of 11 µs for each scan over 4.25 mm in depth. After applying a medium to the biofilm at room temperature, phase-sensitive Doppler OCT was used to scan the biofilm every 5 minutes for an hour (13 times). The temporal fluctuation of the amplitude and phase of the reconstructed OCT images by the internal dynamics of biofilm was observed. Scattered light from a biofilm forms a speckle pattern as the size of S. aureus is smaller than the wavelength. The biofilms were not completely dispersed in our measurement time window after submerging the biofilms in a medium. The interaction between the random walkers and the biofilm demonstrated the coexistence of two distinct statistics (i.e., random walk in free space and biofilm).

Results
We investigated the phase variances of OCT images acquired from biofilms after being submerged in a medium. Here, we assumed that the development of the dynamics in the biofilm was thermally driven, as S. aureus lack flagella and are not motile. The structural deformation induced by applying medium is shown in Fig. 1(b) and (d) and the dynamic characteristics were first investigated by the phase variance. Second, the anomalous diffusion characteristics were further investigated to demonstrate the correlation between the phase variance and the anomalous diffusion characteristics of biofilms. The slow volumetric expansion of the biofilm after being submerged in a medium, and the development of dynamics, were observed by the phase-sensitive Doppler OCT system. Temporally acquired images were coherently reconstructed and the statistics of the phase angles of complex-valued images were acquired. The images of the phase-sensitive Doppler OCT were linearly decomposed into static and dynamic components and the anomalous diffusion characteristics of the corresponding components were investigated. The relationship between the phase variance and the anomalous diffusion characteristics will be discussed and a numerical model to describe the ballistic dynamics of the biofilm will be introduced that can potentially determine whether the biofilm under intervention is undergoing dispersal or not.

Phase-sensitive dynamic measurements from biofilm
Coherently reconstructed images acquired by phase-sensitive Doppler OCT carry dynamic information. The conventional parameters of Doppler OCT to estimate the dynamics of a target are the time-averaged and the variance of phase displacements [66]. Acquired timeaveraged phase displacements and phase variances of each voxel over 100 frames from the phase-sensitive OCT image were converted into the time-averaged Doppler velocity (v z ) and the motility (root-mean-square of the Doppler velocity of a voxel over 100 frames, v rms ) by Eq. (1), respectively. The time-averaged Doppler velocity and motility distributions of biofilm shown in Fig. 2(a) can be induced by various mechanical processes such as biofilm expansion and dispersal [16,20,21,24,67]. The initial osmotic pressure within the biofilm induced a systematic swelling motion which exhibited an average Doppler velocity during the measurement time (∼2 s) of 0.75 µm/s. After an hour, the distribution of the time-averaged Doppler velocity had a mean value of 0, which indicated that the swelling motion was minimized. In some extreme cases like in Fig. 2(a), the biofilm structure had cracks from structural deformation. The time-averaged Doppler velocity and motility distributions in Fig. 2(b) can be fitted by a linear combination of two distinct distributions that suggest the system has two dynamic components with distinct mean v rms values. As the measurements were performed on Gaussian background noise in Eq. (5), any dynamic components of v rms close to the noise level would not be measurable. We assumed the static components of the biofilm are in this regime. The motility distribution in Fig. 2(b) was approximated by a two-term Gaussian distribution whose local maxima of likelihood are defined as v S and v D , which are the means of the static and dynamic components, respectively. The v rms distributions were measured on the Gaussian noise floor, so the v rms distribution of the static components is initially close to the Gaussian distribution with weakly developed dynamics. The PDF of the motility can be expressed as where p S and p D are fractional portions of static and dynamic components (p S + p D = 1). The fractions were acquired by dividing the volume of static and dynamic components by the entire biofilm shown in Fig. 3(a). The volume estimations were performed by counting voxels. Two components were numerically separated by the threshold v Th , which is defined as an interception point between the motility distributions of t = 0 and t = 60 min, and σ S and σ D , which are the widths of each Gaussian distribution. The local maximum of the dynamic component v D is not always well-defined when the PDF shape changes minimally between the initial and final measurements. However, the initial PDF has a well-defined Gaussian-like shape with a time-independent mean and width, so the dynamic component can be expressed as an arbitrary function. In this case, the average motility of biofilm is dominated by the static component, and the contribution of the dynamic component is minimal. The threshold defined in Fig. 2 is independent of the fraction parameters p S and p D . The interception of the motility distribution in Fig. 2 is at p t=0 (v Th ) = p t≠0 (v Th ) and the v th values did not alter much during the transition process, so we assumed the v th is a constant. The motility distribution in Fig. 2(b) was acquired for 50 independent biofilms and the corresponding v th values were estimated. With this condition, the relationship can be derived from Eq. (7) by the by assuming p S ∼ 1 and p D ∼ 0 when t = 0, and v S and σ S are time-independent. The v Th value is independent of the fractional coefficients and only depends on the dynamic parameters of the static and dynamic components. Theoretically, the initial measurement should have p D = 0 at t = 0, however, practical measurements take time to commence data acquisitions, and the internal dynamics were weakly developed during the moment.

Relationship of biofilm dynamics and super-diffusive characteristics
From Fig. 2(b) and Fig. 3(a), two distinct dynamic components were observed in biofilms. From this information, we acquired the PDFs of corresponding components shown in Fig. 3(b) and derived the correlation between the motility and the anomalous diffusion characteristics. The estimation of the driving force causing biofilm dynamics is challenging, but the biofilm volume expansion might be one of the consequences of the driving force as the biofilm volume expansion and the average motility show a linear correlation in Fig. 3(c). Two-term Gaussian components in Eq. (4) imply that the average motility is determined by the coefficients of each component. The relationship between the average motility and the fractional portion was examined by the average motility obtained from Eq.
The fractional volumes of each component are denoted as V S and V D for the static and dynamic components, respectively. The volumes of the static and dynamic components varied over time. The differences of the portions of each component at the time point t = 0 and t = 60 min can be estimated by defining where V tot is the total volume of biofilm, and subscripts i and f represent initial (t = 0 min) and final (t = 60 min) data points, respectively. By substituting Eq. (10a-d) into Eq. (9) with the constraint condition V tot = V S + V D , the fractional portion change (∆p D and ∆p S ) and the average motility change between initial and final measurements (v rms,f − v rms,i ) have the following linear relationship by assuming the initial fractional portion of the dynamic component is small (p D ≈ 0 and p S ≈ 1 at t = 0) where V D,f and V S,f are the fractional volumes of the dynamic and static components at t = 60 min, respectively. Equation (11) and the measurements from 50 independent biofilms shown in Fig. 3(d) show the biofilm motility is determined by the proportion of the static and dynamic components, and that the mean motilities of the static and dynamic components are reproducible. Yet, the estimation of the dynamics by the mean (v z ) and variance of Doppler velocity (v rms ) in Fig. 2(b) are based on the assumption that the statistics of the momentary displacements would be a normal distribution, however, the momentary Doppler velocity distribution v z (not a time-averaged value) showed anomalous distributions in Fig. 3(b). To investigate the anomalous diffusion characteristics, the PDFs of the Doppler velocities of the static and dynamic components were acquired. The biofilm geometry and the corresponding dynamics after submerging the biofilm in a medium are shown in Fig. 3(a). The development of the biofilm dynamics induced the enhancements of the power-law tail in the PDF. The increasing power-law tail of the entire biofilm indicates that the biofilm dynamics become more super-diffusive after the biofilm is submerged in a medium. The power-law tails of the dynamic components are close to the Cauchy diffusion (α ∼ 1, ballistic), while the static components show more Gaussian-like diffusion (α ∼ 2, diffusive). The average motility increment and the Lévy α values in Fig. 3(e) and (f) show a linear correlation, while the dynamic component of higher motility is closer to the Cauchy distribution. The Lévy α values of the static components in Fig. 3(f) do not show noticeable changes according to the motility change.

Relationships between α, γ, and CGMM approximation
The morphological transition of biofilm and the corresponding dynamics were investigated. The interpretation of the dynamics using the associated motility (phase variance) simply represents linear decrement and increment of the static and dynamic components by Eq. (11). However, the dynamic picture of anomalous diffusions in Fig. 3(f) implies that v th only works robustly for the two-term Gaussian approximation in Eq. (7), as the power-law shape of the dynamic components varies. The subtle changes in the power-law tails observed in the PDFs would require more than strictly a linear relationship as the likelihood of the power-law distribution is varied in logarithmic scales. To analyze and estimate the amount of ballistic components, the linear combination of diffusive and ballistic distributions are introduced and the corresponding power-law shapes are analyzed. The result shown in Fig. 3 is consistently indicating the expansion of biofilm, the development of biofilm dynamics, and that more ballistic diffusions are correlated when assuming biofilms had not significantly grown during the relatively short measurement time. During the process, the power-law tails of the PDFs were altered according to the increment of the motility. The interpretation is that the alternations of a Lévy distribution have occurred by nonlinear and fractional-dimension processes, and that the likelihood of the power-law tail component of PDF increases exponentially according to the Lévy α parameter (∼∆ϕ −α−1 ) [29,47]. The analytic derivation of the relationship between the variables demonstrated in Fig. 3 is challenging, so the distributions were approximated by the Cauchy and Gaussian mixture model (CGMM, Fig. 4(a)) [68]. CGMM is the linear combination of the Cauchy and Gaussian distributions with the same scaling parameter γ. Note that the power-law tail shape is sensitive to the linear coefficients of CGMM rather than the scaling parameter γ. Here the Cauchy and Gaussian distributions were acquired with the same scaling parameter γ from Eq. (3) by setting α = 1 and 2, respectively. The equation of the PDF (p(∆ϕ)) by CGMM can be written as where the coefficient has a value 0 ≤ c ≤ 1. The coefficient (1 -c) and c represent the portion of ballistic and diffusive components, respectively. By varying the c value, Lévy α and γ values of the distribution were acquired. From the α and γ table of CGMM, the Cauchy and Gaussian coefficients of the experimental data were estimated. The coexistence of rare ballistic walks and local diffusion is the unique feature of the Lévy walk, and this can be approximated by the probabilistic linear combinations of the Cauchy and Gaussian distributions. The method can simplify the interpretation of Lévy-like distributions we observed from the internal dynamics of biofilm. The shape of CGMM is similar to the PDFs shown in Fig. 4(a), as the observation is statistics of photon ensembles backscattered from The p E and p T parameter distributions acquired by fitting experimental data from 50 independent biofilms using Eq. (11), and the functional shape classifications. The grey region (near 0) represents the quasi-equilibrium state, the red region represents the rapid dispersal process, and the green region represents a less-susceptible process. When p E = p T , the process shows minimal alternation of the E 0 and T 0 portions. ballistic and diffusive dynamic particles, whose displacement distributions are assumed to be Cauchy and Gaussian distribution, respectively. Furthermore, the CGMM may represent the particle dynamics of the trapped phase or escaped phase while interacting within the biofilm ( Fig. 1(a)). All permutations of linear coefficients and widths of the CGMM were numerically generated and the optimized α and γ values were acquired ( Fig. 4(a)). The PDFs acquired from the experiment were converted into α and γ. To characterize the temporal evolution of power-law tail shown in Fig. 4(b) and (c), two coupled equations are proposed as follows The coupled equations describe the number of dynamic particles interacting within the biofilm. Here, T t represents the number of particles trapped inside of biofilm and E t represents the number of escaped particles at time t. The probabilities p E and p T are escaping and trapping probabilities while particles are colliding with the biofilm surface. The packed structure of biofilm initially releases dynamic particles when the biofilm is immersed in a medium. Also, the escaped particles can be re-trapped while colliding with the biofilm. The probabilities of escaping and trapping are assumed to be time-independent, and only the portion of trapped and escaped particles varies temporally.
The terms T t and E t can be decoupled by the matrix diagonalization and the relationships become The biofilm reaches equilibrium when p T E 0 = p E T 0 . The total number of the particles is conserved by the relation T t + E t = T 0 + E 0 . The dynamic ranges of T t and E t are determined by p T and p E .
CGMM could simplify the interpretation of the developments of the power-law tail and be used to derive the integrated dynamic characteristics. The temporally increasing biofilm motility altered α and γ values in Eq. (3). The corresponding coefficients T t /(T t + E t ) and E t /(T t + E t ) of the CGMM were numerically acquired, and the coefficients were fitted to the model function in Eq. (14) by the minimum of residual sum of squares (RSS) in Fig. 4(d) and (e). The p E and p T values were acquired by fitting T t and E t to the model function in Eq. (6). The curvature of the exponential decay was determined by the exponent term 1 -p E -p T . The average value of the initial coefficient 100 × T 0 /(T 0 + E 0 ) was 58.3 ± 6.2, which represents that the portion of the trapped and escaped particles of all biofilms was similar at the initial measurements. The different exponential shapes shown in Fig. 4(e) are classified according to p E and p T and further demonstrated in the colored regions in Fig. 4(f). The biofilms in the grey region showed a quasi-equilibrium state by the slow exponential decay (small p E and p T values). In contrast, the biofilms in the red region, where p T < p E and p E is larger than 1%, showed a rapid dispersal process with a large dynamic range that is physically consistent, thus the dynamics are dominated by escaping particles. When p E = p T , the dynamics are close to the static state as the process is occurring around the equilibrium point p T E 0 = p E T 0 , and the speed of convergence is determined by p E and p E . The biofilms in the green region showed that they were less susceptible to this stimulation induced after being immersed in the medium, thus the escaping rate was smaller than the trapping rate, and instantly released particles were re-absorbed into the biofilm. The most rapid dispersal and least susceptible biofilm dynamics were observed when an antibiotic was applied with the media, which might describe biofilm breakdown and reduced antibiotic susceptibility [7,69].

Discussion
Biofilm formation is a common strategy of bacteria for environmental adaptation. Understanding and controlling biofilm formation and dispersal mechanisms would benefit many industrial and medical applications. Dynamic characteristics of biofilms are still an under-explored area even though various biofilm-forming strains have common dynamic phases such as adhesion, biofilm formation, and dispersal. Each dynamic phase has complicated nonlinear biochemical interactions and the characteristics of the associated dynamics have not been fully understood. Investigating dynamic biofilm properties has faced many technical challenges, one of which is acquiring statistically robust data from diverse and heterogeneous biofilm structures. Combining OCT imaging and phase-sensitive Doppler detection of dynamic light scattering can generate wide-area and label-free motility maps derived from the complex-valued OCT images. The PDFs of phase displacements acquired from the OCT images show a power-law tail distribution which implies a super-diffusion process. Biofilms with higher motility had more super-diffusive PDFs and more volumetric expansion, which is consistent with the biofilm dispersal dynamics. The conventional parameter that has been used to measure whether a biofilm is in a state of dispersal has been to measure the biofilm thickness over a longer time period, however, the technique presented in this study observes the anomalous diffusion of a biofilm, which can be used to determine the response almost instantaneously. For this study, observations required relatively thick biofilms to collect robust statistics from more voxels. Thinner biofilms could still be interrogated with this approach, but would require longer sampling time or a higher spectral resolution from the spectrometer in order to perform the analysis with a smaller voxel size. Finally, following the application of media to the biofilm, measurements were conducted on the dish exposed to air, so the time window for measurements was limited due to evaporation over time. Performing these dynamic biofilm measurements within enclosed cells would provide a more controlled and stable environment for further studies to observe a strong correlation between biological responses of chemical components in a medium and the dynamics of biofilm driven by quorum sensing.
Lévy distributions were introduced to elucidate the dynamic information from the PDFs of phase displacements acquired by phase-sensitive Doppler OCT. A Lévy distribution describes an anomalous diffusion process. Interpretations of the alternations of the power-law tails are related to fundamental questions, but also provide insights into the dynamic responses of systems [52,53]. To derive a simple dynamic model for morphological transitions that occur in a biofilm, a non-motile bacterial species was introduced. S. aureus is a bacterial species commonly found on the skin, and opportunistically causes chronic infections in elderly or immunocompromised patients by forming biofilms on the surfaces of medical implants [70][71][72]. Biofilm formation is problematic as the infection of such patient groups can be fatal, so various medical precautions and methods for detecting and preventing biofilm formation have been developed. S. aureus biofilms are mechanically stable and the motility that develops after submersion in various reagents showed dynamic changes in the structural heterogeneity. The static and dynamic components could be separated by a threshold value defined as an intersection in Fig. 2(b). The threshold value works robustly with the two-term Gaussian approximation and the separated components show the motility contrast in Fig. 3(a). The PDFs of the Doppler velocities collected from the static and dynamic components showed distinct PDF shapes. The temporal PDF shapes of the Doppler velocities show increasing power-law tails. The development of the power-law tail reflects that the system is more ballistic and less diffusive. To interpret the temporal power-law tail alternations, the Lévy α and γ values were converted into the linear coefficients of the CGMM. The temporal evolution of the portions of trapped and escaped particles provides an intuitive physical picture and enables the estimation of the interaction parameters between dynamic particles and the biofilm. This method provides dynamic biophysical parameters that may be used to assess antibiotic susceptibility for medical applications. Additional analytical characteristics of anomalous distributions can be further explored. For instance, the parameters of CGMM and the Lévy distribution were numerically derived, however, analytic derivations could further validate the model and demonstrate the physical significances of the parameters. In addition, the two competing factors introduced in Eq. (14) can be interpreted as biofilm growth and dispersal factors that can determine biofilm phenotypes. CGMM has been used to approximate anomalous distributions, as analytic approaches are challenging due to the fractional operators. While the CGMM approximation can be validated to perform robustly, the analytic investigations of the anomalous distributions would be more efficient. Also, detecting the power-law PDF and fitting with CGMM can estimate the portion of decorrelated photons of the OCT system that can determine the stability of OCT system performances [73][74][75].

Conclusion
This study presents a method to measure the dynamic responses of biofilm to environmental perturbations using images acquired by phase-sensitive Doppler OCT, and by numerical estimations of the dynamic characteristics. The method demonstrates not only the structural deformations of biofilm, but also the corresponding development of anomalous diffusion characteristics. Such dynamic characteristics were numerically estimated by MLE of the Lévy distribution and the experimentally acquired PDFs. For further analysis, experimentally acquired PDFs were approximated by the linear combinations of Cauchy and Gaussian distributions and the development of the power-law tails were investigated numerically. A sessile bacteria model was used to reject active transfer dynamics such as swarming. This technique can be a method to further investigate biofilms with various combinations of bacteria species exhibiting symbiotic or competitive relationships, to study the dormant and active phases for motile bacteria, and to investigate bacterial biofilm responses to various forms of antibiotic or electromagnetic energy (light, plasma) treatments.