Time delay and velocity estimation methods have been widely studied subjects in the context of signal processing, with applications in many different fields of physics. The velocity of waves or coherent fluctuation structures is commonly estimated as the distance between two measurement points divided by the time lag that maximizes the cross correlation function between the measured signals, but this is demonstrated to result in erroneous estimates for two spatial dimensions. We present an improved method to accurately estimate both components of the velocity vector, relying on three non-aligned measurement points. We introduce a stochastic process describing the fluctuations as a superposition of uncorrelated pulses moving in two dimensions. Using this model, we show that the three-point velocity estimation method, using time delays calculated through cross correlations, yields the exact velocity components when all pulses have the same velocity. The two- and three-point methods are tested on synthetic data generated from realizations of such processes for which the underlying velocity components are known. The results reveal the superiority of the three-point technique. Finally, we demonstrate the applicability of the velocity estimation on gas puff imaging data of strongly intermittent plasma fluctuations due to the radial motion of coherent, blob-like structures at the boundary of the Alcator C-Mod tokamak.

## I. INTRODUCTION

Velocity estimation plays a crucial role in a wide range of scientific and technological fields. Its applications extend from analyzing turbulent systems, such as atmospheric phenomena^{1} and magnetized plasmas, to clinical ultrasound scanning,^{2} and to exploring disciplines such as astrophysics and space sciences.^{3} Moreover, it plays a key role in enhancing the effectiveness of radar systems and improving communication technologies.^{4} Many different techniques have been studied, conditioned by the field of application and the experimental setup. In this contribution, we are addressing the problem of estimating the velocity of localized structures in an intensity field recorded in a two-dimensional plane.

One particular example of such a system is fluctuations due to field-aligned filaments at the boundary of magnetically confined plasmas. These are high-pressure structures that are localized in the plane perpendicular to the magnetic field, thus commonly referred to as blobs. In many experimental setups, a gas-puff imaging diagnostic measures light emission from these coherent structures, which move radially out of the plasma column.^{5–7} Different approaches have been employed for velocity estimation in this setup. Most commonly, the velocity is inferred from time-delay estimation (TDE) from the measured signals. TDE methods range from cross correlation techniques^{8–18} to wavelet methods^{19} and dynamic programming.^{20–24} Alternative approaches to velocity estimation, not relying on cross correlations, encompass Fourier analysis,^{25–29} optical flow continuity,^{30,31} blob tracking,^{32–35} and machine learning based algorithms.^{36–39}

In its simplest form, the standard method estimates the velocity of propagation between two measurement points based on the distance between those points and the time delay between the signals recorded. In fluctuating media, this time delay is commonly estimated as the time lag that maximizes the cross correlation function between the signals. In the following, we will refer to this as the two-point method. This method is inaccurate if the velocity of propagation is not aligned with the measurement points or when the propagating structures are tilted.^{28,29,40} The commonly used Fourier analysis methods are subject to the same limitations as this two-point estimation method.^{14,26–28} Improved methods for velocity estimation have been developed using spatial cross correlation functions, but these are limited to high image resolution diagnostics.^{17,30,41–43}

The key idea underlying the method lies in the realization that the estimated time delay between two measurement points is not determined solely by the component of the velocity vector in the direction of those two points. Instead, it depends on both components of the velocity vector. This phenomenon was previously observed, particularly in Ref. 28, where synthetic data were used to estimate velocities using a two-point approach. It was found that the input velocities were off by a factor that could be corrected using the perpendicular component of the velocity; however, this observation was not further developed into a velocity estimation method. In Ref. 29, similar relations were obtained in the case of propagating wave-like structures and were applied to develop a velocity estimation method based on cross-phase analysis. In this work, we rigorously demonstrate that the same principles can be applied to the propagation of localized pulses and that the results are valid for any distribution of pulse amplitudes and density of pulses. In addition, we benchmark the method against synthetic data to confirm that it accurately determines velocities for any direction and to explore its applicability dependence on spatial resolution.

In Sec. II, we present an improved three-point method to estimate the velocity of fluctuation structures moving in a two-dimensional plane from simple geometric considerations. In Sec. III, a stochastic model describing the fluctuations as a superposition of Gaussian pulses moving with constant velocities in the two-dimensional plane is introduced. This demonstrates that the three-point method can be used with the time delays estimated from cross correlation functions and that it is independent of the distribution of pulse amplitudes as well as the density of pulses. This is verified by the analysis of synthetic data from realizations of this process in Sec. IV. In Sec. V, the three-point method is applied to measurement data from gas puff imaging experiments at the boundary of the Alcator C-Mod tokamak. A discussion of the TDE method and conclusions of these investigations are presented in Sec. VI. Finally, Python implementations of both the velocity estimation methods and the synthetic data generation described in this paper are openly available on GitHub.^{44,45}

## II. TIME DELAY ESTIMATION

Consider a front propagating along its normal in a two-dimensional plane with horizontal velocity component *v* along the *x*-axis and vertical velocity component *w* along the *y*-axis. The perturbation is measured at three spatially separated points *P*_{0}, *P*_{x}, and *P*_{y} as illustrated in Fig. 1. In this arrangement, *P*_{x} is horizontally separated from *P*_{0} by a distance △_{x}, while *P*_{y} is vertically separated from *P*_{0} by a distance △_{y}. We denote the velocity vector angle from the horizontal axis by *α*.

*v*and

*w*from the signals measured at

*P*

_{0},

*P*

_{x}, and

*P*

_{y}. Consider first the recordings at

*P*

_{0}and

*P*

_{x}. Let

*τ*

_{x}be the time delay between the measurements of the front at these positions. In the standard two-point TDE method, this time lag is interpreted as the time taken for the pulse to traverse horizontally the distance △

_{x}between points

*P*

_{0}and

*P*

_{x}, that is,

*τ*

_{x}= △

_{x}/

*v*. From this, the two-point estimate for the horizontal velocity component

*v*is given by

*τ*

_{y}between the signals measured at

*P*

_{0}and

*P*

_{y}, the vertical velocity component is estimated as

However, this interpretation of the time lags *τ*_{x} and *τ*_{y} is obviously flawed.^{28,29,40} To illustrate this, consider the special case of a horizontal front that moves strictly vertically, rendering *v* = 0. Such a pulse passes through the measurement points *P*_{0} and *P*_{x} simultaneously, resulting in a time lag *τ*_{x} = 0. This would then lead to an infinitely large estimated horizontal velocity component $v\u03022$. More generally, if the pulse moves at a slight angle to any coordinate axis, the two-point method will give a severe overestimate of the corresponding velocity component.

*P*

_{0},

*P*

_{x}, and

*P*

_{y}, it is straightforward to give correct estimates of the two velocity components. As the pulse moves, it will first be measured at

*P*

_{0}and some time later at

*P*

_{x}, as seen in Fig. 1. The structure will have traveled a distance $\u25b3xcos\alpha =v\u25b3x/(v2+w2)1/2$ along its direction of motion between its arrivals at

*P*

_{0}and

*P*

_{x}. The time lag

*τ*

_{x}corresponds then to the time △

_{x}cos

*α*/

*u*it takes for the pulse to travel this distance,

*P*

_{0}is required. A similar geometrical consideration gives the distance the structure travels as $\u25b3ysin\alpha =w\u25b3y/(v2+w2)1/2$, where △

_{y}is the vertical separation of the measurement points. The resulting lag time is then given by △

_{y}sin

*α*/

*u*, or in terms of the velocity components,

*v*and

*w*, leading to the velocity component estimates,

*P*

_{0}and

*P*

_{x}, the time lag

*τ*

_{x}vanishes and so does the estimated horizontal velocity component from Eq. (3a). A similar three-point method has previously been discussed, recommending the use of the cross-phase from the cross-power spectral density to estimate the time delays.

^{29}

These considerations demonstrate that the two velocity components of a sharp front or a plane wave can be accurately estimated by using three non-aligned measurement points. It is straightforward to generalize this to a localized pulse that is symmetric with respect to its direction of motion, which gives the same expressions for the velocity components. However, in this case, the distance between the measurement points must be smaller than the structure size. Moreover, in turbulent flows, there can be an overlap of structures and a distribution of amplitudes, sizes, and velocities. This motivates a statistical treatment, which is presented Sec. III.

## III. STOCHASTIC MODELLING

It is generally anticipated that the TDE presented in Sec. II for the case of a single front can be generalized to a fluctuating state with multiple structures. For a statistically stationary process, the time lags *τ*_{x} and *τ*_{y} are then typically taken as the correlation times that maximize the cross correlation functions for measurements at the positions *P*_{0}, *P*_{x}, and *P*_{y}. In this section, it is demonstrated that Eq. (3), indeed, gives exact results in the case of a superposition of uncorrelated Gaussian pulses having all the same size and velocity components. The derivation provided here will be schematic, with a more detailed and comprehensive study presented in a separate publication.

*K*uncorrelated pulses,

*a*

_{k}arrives at

*x*= 0 and

*y*=

*y*

_{k}at time

*t*=

*t*

_{k}, and where

*y*

_{k}and

*t*

_{k}are uniformly distributed random variables. The process has duration

*T*, and the vertical domain size is

*L*. The pulses move with velocity components (

*v*,

*w*) and have horizontal and vertical sizes

*ℓ*

_{x}and

*ℓ*

_{y}, respectively. All pulses are taken to have the same functional form

*φ*, which is here assumed to be Gaussian,

^{46–51}

*a*⟩ is the average pulse amplitude, and

*τ*

_{w}=

*T*/⟨

*K*⟩ is the average waiting time between pulses. The factor

*ℓ*

_{x}/

*vτ*

_{w}is the ratio of the radial transit time

*ℓ*

_{x}/

*v*to the average waiting time and, therefore, determines the degree of temporal pulse overlap. Similarly,

*ℓ*

_{y}/

*L*determines the degree of spatial pulse overlap in the vertical direction. It can furthermore be shown that the variance of the process is given by

_{rms}/⟨Φ⟩ becomes large when there is little overlap of pulses.

^{46–51}

_{x}, △

_{y}, and △

_{t}are lags in space and time. Performing this average using the Gaussian pulse shape leads to the following cross correlation function:

_{t}maximizing

*R*

_{Φ}in Eq. (10) for fixed spatial lags △

_{x}and △

_{y}is given by

*τ*

_{x}and

*τ*

_{y}maximizing the cross correlation functions at the fixed spatial lags of the given setup can be estimated. In particular, for △

_{y}= 0, the cross correlation function is maximum for time lag

*τ*

_{x}=

*v*△

_{x}/(

*v*

^{2}+

*w*

^{2}), while for △

_{x}= 0, the cross correlation function is maximum for time lag

*τ*

_{y}=

*w*△

_{y}/(

*v*

^{2}+

*w*

^{2}). This is identical to the results in Eq. (3) obtained by geometrical considerations, thus leading to the same expressions for the velocity components. It is also straightforward to obtain velocity estimates from any two time delay estimates, meaning that the three points used by the method need not be orthogonal or aligned with the coordinate grid. Finally, it is straightforward to show that the time of maximum cross correlation given by Eq. (12) agrees with the cross-phase TDE suggested in Ref. 29.

The results from this stochastic modeling cannot be overemphasized. First, it shows that the three-point calculation of pulse velocities obtained from simple geometric considerations is exact in the case of a superposition of uncorrelated Gaussian pulses which all have the same size and velocity. Second, this estimate is independent of the distribution of pulse amplitudes as well as the degree of pulse overlap or the density of pulses. Finally, it confirms that the cross correlation function can, indeed, be used for TDE.

## IV. NUMERICAL SIMULATIONS

In order to confirm the theoretical predictions above, we performed numerical realizations of the stochastic process. The Gaussian pulses are taken to be symmetric, *ℓ*_{x} = *ℓ*_{y} = *ℓ*, and the total speed $u=(v2+w2)1/2$ is fixed while the ratio of the velocity components *v*/*w* varies between realizations. The vertical domain size is set to *L* = 10*ℓ*, and periodic vertical boundary conditions are implemented. The amplitude is, for simplicity, taken to be the same for all pulses, and the total duration of the process is *uT*/*ℓ* = 10^{3} with *K* = 10^{3} pulses in total. The resulting fluctuations are measured at a fixed reference point *P*_{0}, and two auxiliary points *P*_{x} and *P*_{y} separated by a pulse size *ℓ* in the *x*-direction and in the *y*-direction, respectively. The sampling time is 10^{−2}*ℓ*/*u*, thus representing typical experimental conditions with high temporal sampling rates but relatively coarse-grained spatial resolution of the measurements, as detailed in Sec. V.

Simulations are carried out for various combinations of input velocities *v* and *w*. For each case, the estimated velocities $v\u0302$ and $w\u0302$ are computed using either the two- or the three-point method as described by Eqs. (1) and (3), respectively. In both cases, the time lags are obtained from the cross correlation functions. The results of the velocity estimation are presented in Fig. 2, where the estimates resulting from the two- and three-point methods are shown with crosses and filled circles, respectively.

As suggested by the stochastic model in Sec. III, the three-point method gives correct estimates for both velocity components in all cases, despite a likely marginal spatial resolution for the measurement points, which are separated by a pulse size *ℓ* in the numerical simulations. However, the two-point method fails as expected from the discussion in Sec. III. When the pulse velocity vector is nearly parallel to one of the coordinate axes, the error for the opposite velocity component can be arbitrarily large. Indeed, the two-point estimate of the velocity vector can, according to Eq. (2), be written as $(v\u03022/u,w\u03022/u)=(u/v,u/w)$, which is represented by the dashed lines in Fig. 2. This is an excellent description of the simulation results and demonstrates the failure of the two-point method when one of the velocity components is much smaller than the other.

In order to quantify this, consider the case where we want to ensure that the relative error from the two-point estimate of the vertical velocity component $w\u03022$ remains below a value *p*. According to Eq. (2b), this means that $w\u03022/w=1+v2/w2<1+p$. It follows that the ratio of the horizontal to vertical velocity components is given by *v*^{2}/*w*^{2} < *p*, that is, when the error on the estimate of the vertical component is small, *p* ≪ 1, the horizontal velocity component is significantly smaller than the vertical component. From Eq. (2a), it furthermore follows that the relative error for the horizontal velocity component is $v\u0302/v=1+w2/v2=1+1/p$. So, for a 10% relative error on the vertical component, *p* = 10^{−1}, the error for the horizontal velocity component is an order of magnitude larger, $v\u03022/v=11$, when using the two-point estimates. Of course, in practical applications, the true velocity components are unknown and there is no way to infer the relative error from the two-point method.

*v*=

*w*and varying spatial sampling distances △

_{x}= △

_{y}= △. Moreover, additive temporally and spatially uncorrelated white noise has been added to simulate measurement noise.

^{52}The results of the three-point velocity estimation method are presented in Fig. 3, where the squared error is defined by

*E*< 0.1) for spatial resolutions in the range 10

^{−2}< △/

*ℓ*< 5. When △/

*ℓ*< 5 × 10

^{−2}, the time maximizing the cross correlation between the signals at different points is of the order or the sampling time d

*t*= 10

^{−2}

*ℓ*/

*u*, which leads to errors in its estimation even though the cross correlation function is interpolated in time. A possible way to solve this is to use more distant measurement points if those exist. On the other hand, when △/

*ℓ*> 5, a pulse detected at one point might not be detected at the other measurement points, meaning that the signals at different points become uncorrelated and it is not possible to accurately estimate the cross correlation maximum, and hence, we obtain the order unity error shown in Fig. 3 for these cases.

A comprehensive numerical simulation study has been performed, testing the three-point method for velocity estimation for various pulse functions and distributions of pulse amplitudes, sizes, and velocities as well as sampling rates and density of pulses. It is generally found to give accurate estimates of the average velocity components, while the two-point method consistently fails to reliably estimate both components. The details of this study will be presented in a separate report.

## V. EXPERIMENTAL MEASUREMENTS

Alcator C-Mod is a compact, high-field tokamak with major radius *R*_{0} = 0.68 m and minor radius *r*_{0} = 0.21 m.^{53–55} The device is equipped with a gas puff imaging diagnostic, which consists of a 9 × 10 array of in-vessel optical fibers with toroidally viewing, horizontal lines of sight, as shown in Fig. 4.^{7,26} The plasma emission collected in the views is filtered for He I line emission (587 nm) that is locally enhanced in the object plane by an extended He gas puff from a nearby nozzle. Because the helium neutral density changes relatively slowly in space and time, rapid fluctuations in He I emission are caused by local electron density and temperature fluctuations. The GPI intensity signals are, therefore, taken as a proxy for the plasma density.^{5,7}

The optical fibers are coupled to high-sensitivity avalanche photodiodes, and the signals are digitized at a rate of 2 × 10^{6} f/s. The viewing area covers the major radius from 88.00 to 91.08 cm and vertical coordinate from −4.51 to −1.08 cm with an in-focus spot size of 3.8 mm for each of the 90 individual channels. The measurements presented here were done during the last operational year of Alcator C-Mod. Numerous diodes were broken, leading to a lack of data for several channels in the 9 × 10 array of measurement points.

The experiment analyzed here was a deuterium fueled plasma in a lower single null divertor configuration with only Ohmic heating. We present an analysis of the plasma fluctuations recorded using the gas puff imaging system for discharge number 1 160 629 026 in the time frame from 1.22 to 1.50 s. This discharge had a plasma current of *I*_{p} = 1.1 MA, an axial magnetic field of *B*_{0} = 5.4 T, and a Greenwald fraction of line-averaged density $n\u0304e/nG=0.27$, where the Greenwald density is given by $nG=(Ip/\pi r02)1020m\u22123$ with *I*_{p} in units of MA and *r*_{0} in units of meters.^{56} Analysis of some of these measurement data has previously been reported in Refs. 57 and 58 with a focus on the fluctuation statistics in the far scrape-off layer.

Cross-field particle transport at the boundary of magnetically confined plasmas is generally attributed to localized blob-like filaments moving radially outward toward the main chamber wall, resulting in intermittent and large-amplitude fluctuations in the plasma parameters.^{57–66} Such radial motion is evident from the raw measurement data time series at different radial positions, presented in Fig. 5. Here, we show the measured line intensity signals at *Z* = −2.2 cm, where each time series Φ(*t*) is normalized so as to have vanishing mean value and unit standard deviation, $\Phi \u0303=(\Phi \u2212\u27e8\Phi \u27e9)/\Phi rms$. Here, the mean value ⟨Φ⟩ and the standard deviation Φ_{rms} are calculated by a running mean over ∼1 ms in order to remove low-frequency trends in the time series. In Fig. 5, several large-amplitude events propagate radially outward with velocities of the order of 1000 m/s.

The velocity of these fluctuations has been estimated with both the two- and three-point TDE methods described above, using the lag times that maximize the cross correlation function for radial and vertical displacements. For each diode view position, the velocity components are estimated based on correlations with all nearest neighbor diode view positions and averaged radially and vertically, as depicted in Fig. 6. The estimated velocities with these two methods are presented in Fig. 7. The crosses correspond to the location of broken diode view positions. No velocity is assigned if either the cross correlation function is not unimodal or the time of maximum cross correlation is less than the sampling time △_{t} = 5 × 10^{−7} s.

The gray shaded vertical stripe in Fig. 7 between the major radii of 88.8 and 89.2 cm is the location of the last closed magnetic flux surface according to magnetic equilibrium reconstruction. This vertical line separates the confined plasma column to the left and the scrape-off layer to the right of this region. The dotted line at ∼91 cm is the location of limiter structures mapped to the gas puff imaging view position.

The results from the two-point method, shown to the left in Fig. 7, suggest predominantly vertical motion through the edge and scrape-off layer regions. However, as anticipated from the theoretical considerations and numerical simulation studies presented above, the results from the three-point method, shown to the right in Fig. 7, reveal that the motion is mainly radial with velocities up to 750 m/s. The two-point method reverses the dominant velocity components and is obviously wrong, as also indicated by the raw data in Fig. 5, which reveals the radial motion of the pulses. The results from the three-point method presented here qualitatively agree with previous analysis of Phantom camera data, which has a much higher spatial resolution and therefore allows more sophisticated correlation analysis and time delay estimation methods.^{67}

## VI. DISCUSSION AND CONCLUSIONS

In this study, we have proposed and evaluated an improved method for estimating the velocity components of fluctuation structures that move in two spatial dimensions. The method utilizes cross correlation functions obtained from recordings of a scalar intensity field at three non-aligned measurement points, self-consistently taking into account the two-dimensional nature of the problem. This is demonstrated to give an exact expression for the fluctuation velocity components in a stochastic model given by a superposition of uncorrelated Gaussian pulses. This analytical calculation manifests that the cross correlation function can be used for the TDE estimation and that the results hold for any distribution of pulse amplitudes and density of pulses. The accuracy of the method is further demonstrated by application to numerical realizations of the process, supporting its usefulness also in the case of random distributions of pulse parameters.

The results from the improved three-point method presented here are compared to a standard two-point method. The latter essentially describes the problem as one-dimensional and only gives an accurate estimate of the velocity when the two measurement points are perfectly aligned with the direction of motion of the pulses. If the pulses move at a slight angle to any coordinate axis, the two-point method will give a severe overestimate of the corresponding velocity component. Moreover, the two-point method always overestimates both velocity components. This tendency for overestimation has previously been observed in turbulence simulations where the velocity field is known.^{68} Even in that case, only the velocity component in the direction of the two measurement points is correctly estimated, while the perpendicular component would artificially be estimated as infinity. Commonly used Fourier methods with wave number-frequency spectra for one spatial dimension suffer from the same shortcomings as the two-point method.^{14,26,28}

Applying the two velocity estimation methods to experimental measurement data from the avalanche photodiode gas puff imaging diagnostic on Alcator C-Mod reveals a striking difference. The two-point method suggests a dominantly vertical motion of the fluctuation structures, while the three-point method demonstrates that the motion is mainly radial as can also be inferred from the raw data time series. The two-point method gives inherently erroneous results, overestimating the horizontal and vertical velocity components and reversing their ratio. The improved three-point method is, therefore, likely to give results that are consistent with other analysis methods and diagnostics. In future work, this improved velocity estimate method will be used to investigate how fluctuations in the boundary region change with experimental control parameters, improved confinement modes, the role of auxiliary heating, and differences between the far scrape-off layer dominated by the motion of blob-like filaments and the closed field line region where wave dynamics is expected to prevail.

The three-point method presented here has been derived for the case that the pulses are symmetric, *ℓ*_{x} = *ℓ*_{y}, in which case tilting has no effect, or for the case that *ℓ*_{x} ≠ *ℓ*_{y} but the pulses are tilted in such a way that its motion is parallel to one of its axes. In the cases that *ℓ*_{x} ≠ *ℓ*_{y} and the pulses are tilted in any other direction, the presented method will have a biased error. This error is equivalent to the so-called barberpole effect.^{18,69,70} A more technical study of the method, which will be published in a separate contribution, shows that the error made by this simplification will be small if *ℓ*_{x}/*ℓ*_{y} does not deviate much from 1 or if the tilting angle is small. In either case, the error will be much smaller than that made by a two-point approach.

In conclusion, this study presents a valuable approach for estimating velocity components in imaging data or stochastic processes characterized by waves or superposed and localized pulses. By estimating two temporal cross correlation functions in perpendicular directions, the method offers a reliable technique that surpasses the constraints of the standard two-point methods extensively used in previous studies. The improved three-point TDE should, therefore, be added to the list of useful methods to analyze imaging data, in general, and gas puff imaging of magnetically confined plasmas, in particular.^{39} As a final note, it is remarked that the same methodology may be applied with conditional averaging instead of cross correlation estimation for the time delays. This will be particularly useful for investigating flows associated with large-amplitude fluctuations.

## ACKNOWLEDGMENTS

This work was supported by the UiT Aurora Center Program, UiT The Arctic University of Norway (2020). Alcator C-Mod data were generated and maintained under U.S. Department of Energy Award Nos. DE-FC02-99ER54512 and DE-SC0014264. Support for J. L. Terry was provided by the U.S. Department of Energy, Fusion Energy Sciences, Award No. DE-SC0014251. Discussions with S. Ballinger and S. G. Baek are gratefully acknowledged. O.E.G. and A.D.H. acknowledge the generous hospitality of the MIT Plasma Science and Fusion Center where part of this work was conducted.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**J. M. Losada**: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **A. D. Helgeland**: Methodology (equal); Visualization (equal). **J. L. Terry**: Data curation (equal); Resources (equal). **O. E. Garcia**: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.