2.1.

Theoretical Derivation of the Shear Flow MPFRAP Model

The initial concentration distribution of unbleached fluorophores immediately following a photobleaching pulse is given by Brown et al.¹

The bleach intensity can be represented as a 3D Gaussian that is a function of the time average of the intensity at the two-photon focal volume center raised to the $b$’th power, and the $1/e^2$ radial and axial dimensions of the two-photon focal volume, ${\omega }_r$ and ${\omega }_z$, respectively:¹

Assuming a properly overfilled back aperture, the radial and axial dimensions of a two-photon focal volume are defined as ${\omega }_r\equiv 2.6\lambda /(2\pi \mathrm{NA})$ and ${\omega }_z\equiv 8.8n\lambda /[2\pi {(\mathrm{NA})}^2]$, respectively, where $\lambda$ is the wavelength of the excitation laser, $n$ is the index of refraction of the immersion media, and NA is the numerical aperture of the lens.¹⁷

The fluorescence intensity at time $t$ generated by a stationary weak monitoring beam that produces fluorescence through an $m$-photon process is given as¹

The original MPFRAP model assumes that fluorescence recovery is due to diffusion only. The concentration profile, $c(x,y,z,t)$, is determined using the diffusion equation and the resultant expression for fluorescence recovery is then

This model introduces a third fitting parameter, which describes the characteristic recovery time due to flow, ${\tau }_v$. The flow speed in this model can be calculated by the relationship $v={\omega }_r/{\tau }_v$.

The previous model assumes that flow is uniform, but many biological and experimental systems are dominated by shear flow whereby the flow speed in one direction (i.e., the $x$ direction) varies with position in another direction (i.e., the $z$ direction). To make MPFRAP applicable to more in vivo systems, we will therefore incorporate shear flow in the convection element by adding a time-and position-dependent coordinate shift to model shear along the axis perpendicular to flow, before the mathematical convolution of the concentration profile with the excitation laser profile.

To begin the derivation, consider the time-dependent concentration profile of unbleached fluorophores evolving only due to diffusion, as derived in Ref. 1 but expressed here in Cartesian coordinates

To incorporate flow in the presence of shear stress, we apply a time-dependent coordinate shift to this distribution representing flow along the $x$-axis which varies with $z$. In the frame of reference of the observer, we define $x^{\prime }=x+v_{x}t$, $y^{\prime }=y$, and $z^{\prime }=z$. Here $v_x=v_o+\gamma z^{\prime }$ where $\gamma$ is the shear rate and $v_o$ is the flow speed in the x-direction at the center of the focal volume. The time-dependent fluorophore concentration now becomes

By substituting Eq. (10) into Eq. (3), we now have

Since $v_x$ is a function of $z^{\prime }$, we must separate and rearrange the terms in the $x^{\prime }$ and $z^{\prime }$ integrals. We will define $J(t)$ as the product of the three integrals in Eq. (11), and we begin by solving the $y^{\prime }$ integral using the integral identity ${\int }_{-\infty }^{+\infty }{e}^{-{\rho }^2{x}^2}\mathrm{d}x=\sqrt{\pi }/\rho (\rho >0)$. Now, $J(t)$ becomes

To begin integrating the $x^{\prime }$ integral we first define $v(z^{\prime })=v_0+\gamma z^{\prime }$, then define $u^{\prime }(x^{\prime }{z}^{\prime })$ with Eq. (13) and observe that $d{u}^{\prime }/d{x}^{\prime }=1$; this will allow us to make a variable substitution to convert $d{x}^{\prime }$ to $d{u}^{\prime }$ as shown in Eq. (14)

Plugging in the $u^{\prime }(x^{\prime }{z}^{\prime })$ and rearranging the terms such that all solely $z^{\prime }$-dependent terms are in the $z^{\prime }$ integral, $J(t)$ now becomes

Now the $u^{\prime }$ integral is in the same form as the previous $y^{\prime }$ integral, so we can apply the same integral identity to now obtain a simplified version of $J(t)$ with only the $z^{\prime }$ integral remaining

Through rearranging and collecting terms and performing integration with the integral identity ${\int }_{-\infty }^{+\infty }{e}^{-(a{x}^2+bx+c)}\mathrm{d}x=\sqrt{\frac{\pi }{a}}{e}^{(b^2-4ac)4a}$, we obtain a final form of $J(t)$,

We can plug this expression for $J(t)$ back into Eq. (11) for the three integrals, and then define $F_o$, the prebleach fluorescence, using the following equation:

We now introduce our recovery time variables using the definitions: ${\tau }_D={\omega }_r^2/8D$, ${\tau }_v={\omega }_r/v_o$, and ${\tau }_{\gamma }=1/\gamma$. Recalling the fact that MPFRAP is usually performed with two-photon processes, we plug in $m=b=2$ to further simplify, and we also recall that $R={\omega }_z^2/{\omega }_r^2$. We now have the final fluorescence recovery equation that accounts for shear stress in the presence of laminar flow

For simplicity, the exponential term in $J(t)$ from Eq. (16) was named $S_n(t)$. Note that in the case where there is no shear stress (${\tau }_{\gamma }\to \infty$) this model reduces to the diffusion-convection model given by Eq. (5). In the case where there is no flow $({\tau }_v\to \infty )$ and no shear stress (${\tau }_{\gamma }\to \infty$), this model reduces to the diffusion only model, given by Eq. (4).

Figure 1 shows the predicted fluorescence recovery curves using the new shear flow model derived above, with typical values of $\beta =0.6$ and $D=60\mu {\mathrm{m}}^2/\mathrm{s}$. In Fig. 1(a), we observe fluorescence recovery curves in the absence of flow at the center of the focal volume ($v_o=0$) (reminiscent of a “whirlpool” exactly centered on the focal volume), at four different shear rates. Figure 1(b) shows a more realistic scenario with non-zero flow speed at the center of the focal volume. As expected, the presence of shear tends to accelerate recovery of fluorescence.

Fig. 1

Mathematical modeling of MPFRAP with shear flow. Theoretical MPFRAP curves calculated at various stress rates using the shear flow model Eq. (21) where $D=60\mu {\mathrm{m}}^2/\mathrm{s}$, and $v_o=0$ (a) and $v_o=300\mu \mathrm{m}/\mathrm{s}$ (b).

2.2.

Monte Carlo Simulation

To generate artificial data upon which to test our new model, we simulated the motion of bleached molecules under different diffusion/flow/shear conditions. As in previous MPFRAP simulations,^10,17 space was discretized into a regular lattice structure with spacing defined by the expected diffusive properties and diffusion was modeled as a random walk on this 3D lattice.^18,19 The lattice spacing was determined by the 3D diffusion equation $⟨r^2⟩=6Dt$, where the diffusion coefficient, $D$, was chosen a priori and $t$ represents the time step. For this study, $D$ was chosen to be $60\mu {\mathrm{m}}^2/\mathrm{s}$, which is a typical diffusion coefficient for fluorescein isothiocyanate (FITC)-bovine serum albumin (BSA) at room temperature. The bleach depth parameter was set to 0.6.

Previously,^10,17 the time step was set equal to 1/1000 of the typical diffusive recovery time for a system with a diffusion coefficient $D$ and with radial and axial focal volume dimensions ${\omega }_r$ and ${\omega }_z$, thus $t_{\text{step}}={\tau }_D/1000$.^10,17 However, for certain combinations of high flow and large shear rate, the half recovery time of the curve may be smaller than ${\tau }_D$, meaning that the time step would be too large to accurately simulate a fluorescence recovery curve. Instead, we account for all possible recovery mechanisms by calculating the time step as 1/1000 of the expected half recovery time of the curve, calculated using Eq. (20) below. Defining the time step as 1/1000 of the half recovery time due to all possible recovery mechanisms instead of 1/1000 of the half recovery time due to just diffusion^10,17 will ensure a small enough time step to accurately simulate a fluorescence recovery curve at all combinations of flow and shear rate

To start the simulation, this lattice was populated with one candidate molecule at every lattice point, creating a uniform distribution of bleached fluorophores, spanning from $-2{\omega }_r$ to $2{\omega }_r$ in the $x$- and $y$-directions, and $-2{\omega }_z$ to $2{\omega }_z$ in the $z$-direction. In our simulations, the NA was 0.8 and hence ${\omega }_r=0.404\mu \mathrm{m}$ and ${\omega }_z=2.27\mu \mathrm{m}$. Next, a “probability threshold” for each candidate molecule was calculated using Eq. (21) below, which is obtained by substituting Eq. (2) into Eq. (1), setting $b=2$ for a two-photon bleaching process, and $c_o=1$, and represents a distribution of the appropriate shape for a given bleach depth parameter $\beta$. Each candidate molecule was assigned a random number between zero and one, using a uniformly distributed random number generator. The candidate was eliminated if the random number was above the threshold. This process was repeated until there were 20,000 bleached molecules on the lattice. There were no constraints on how many molecules could be at a single node.

Once the initial distribution was created, all molecules were allowed to take one step in a random direction. A single molecule can move in one of six directions; either the positive or negative $x$-direction, positive or negative $y$-direction, or positive or negative $z$-direction. Thus, each molecule was assigned a random integer from 1 to 6, using a uniformly distributed random number generator, where each integer is assigned to one of the six possible directions. The length of each of these steps is defined as $L=\sqrt{6Dt}$. To simulate shear flow, each molecule experienced a horizontal displacement in the $x$-direction after each step, defined as $x_{\text{bias}}=(v_x+\gamma z)t$, where $v_x$ is the flow speed at $z=0$; $\gamma$ is the shear rate; $z$ is the $z$-coordinate of the particle after its random step of length $L$; and $t$ is the time between each step. Once all molecules have taken their random step and have been horizontally displaced, the fluorescence signal from the distribution of molecules at that time step is calculated.

In a physical MPFRAP experiment, the fluorescence of unbleached molecules is monitored; however, simulating an infinitely large volume of molecules is less feasible and more computationally expensive in simulation. For this reason, we are simulating a finite number of bleached fluorophores. Since we are monitoring the bleached fluorophores, we need to calculate the “missing fluorescence” that these fluorophores would have produced had they not been bleached. We can then use the missing fluorescence to determine the fluorescence of the remaining unbleached molecules, which will be used to generate the MPFRAP curve.

Calculating the missing fluorescence of the unbleached molecules can be done by re-expressing the integral in Eq. (3) as a sum of monitor intensity, given by Eq. (2) with $b\to m$, over all bleached fluorophore locations $(x_i,y_i,z_i)$, yielding

Now that we have the missing fluorescence of the bleached fluorophores, we obtain the fluorescence of the unbleached fluorophores, $F(t)$, by recognizing that at any given time after the bleaching event, the fluorescence of the unbleached molecules will be equal to the difference between the prebleach fluorescence ($F_o$) and the missing fluorescence of the bleached molecules, thus $F(t)=F_o-F_{bl}(t)$. Therefore, to produce $F(t)$ from $F_{bl}(t)$, we first need to deduce $F_o$ from only the information carried in the $F_{bl}(t)$ curve. We can deduce $F_o$ specifically from $F_{bl}(0)$ by first taking the $t=0$ limit of Eq. (19),

And substituting this into the $t=0$ expression $F(0)=F_o-F_{bl}(0)$, then solving for $F_o$:

Once the fluorescence recovery curve $F(t)$ is simulated, additional noise is introduced to the curve to represent experimental in vitro and in vivo MPFRAP curves more closely.^1,10,17 Specifically, the poissrnd() function in MATLAB was used to generate a Poisson distributed random number with mean of one, and a distribution width of $\sim 0.03$. This random number was multiplied by each $F(t)$ data point, producing final simulated FRAP curves with a relative noise fraction of 3%, which closely mimics the relative noise found in in vitro experiments (see Fig. 2).

Fig. 2

Example Monte Carlo simulated recovery curve. Fluorescence recovery curve generated with MPFRAP molecular dynamics algorithm in the presence of shear flow. Data are shown with the addition of 3% relative noise, and overlayed is the calculated curve from the mathematical model in Eq. (19).

The simulated fluorescence recovery curves are then fit to the three fitting models (diffusion-only, diffusion-convection, and shear flow) to extract the diffusion coefficient. A typical simulated fluorescence recovery curve is shown in Fig. 2, with the addition of 3% noise and the corresponding theoretical curve using the shear flow model.

2.3.

Fitting Algorithm and Data Analysis

Simulated data were fit to the three fitting models using the Levenberg-Marquardt least squares fitting method implemented in MATLAB (The MathWorks, Natick, Massachusetts, United States) using the lsqcurvefit function. This function requires an initial guess for each fitting parameter as an input. These seed values were calculated as follows.

First, a seed value for the bleach depth parameter $\beta$ was determined by evaluating the shear-flow MPFRAP equation [Eq. (19)] at $t=0$ and comparing it with the first data point of the fluorescence recovery curve after photobleaching, $F(0)$. The seed value of $\beta$ is the value of $\beta$ that satisfies Eq. (25) for a given recovery curve, where $N=10$ for computational efficiency.

Next, we generate seed values for ${\tau }_D$, ${\tau }_v$, and ${\tau }_{\gamma }$. For each term we do this by assuming that the recovery is dominated entirely by that process. Starting with the seed value for ${\tau }_D$, we assume all recovery is due to diffusion. In an automated fashion, we obtain the data point that is closest to or at the half-recovery fluorescence value and define that corresponding time as ${\tau }_H$. The seed value for ${\tau }_D$ is defined as the half-recovery time of the MP-FRAP curve, ${\tau }_H$.

The seed value for ${\tau }_v$ is generated by assuming that recovery is only due to flow. In this case, we take the limits ${\tau }_D\to \infty$ and ${\tau }_{\gamma }\to \infty$ in Eq. (19), which leaves us with an expression for the flow-dominated fluorescence recovery:

We then make the substitution $x^{\prime }=(t/{\tau }_v{)}^2$ and plot $F(x^{\prime })/F_o$ as a function of $x^{\prime }$. The seed value for ${\tau }_v$ is then obtained by finding the half-recovery time of this curve, which is denoted $x_{1/2}^{\prime }$. This value was determined to be $x_{1/2}^{\prime }=0.3625$. Plugging back in for $x^{\prime }$, we can solve for an expression for the seed value of ${\tau }_v$ as follows:

Last, we determine a seed value for ${\tau }_{\gamma }$ by taking the limits ${\tau }_D\to \infty$ and ${\tau }_v\to \infty$ in Eq. (19), which leaves us with an expression for shear flow dominated fluorescence recovery:

Making a similar substitution and using the same mathematical steps we used to calculate ${\tau }_{v_{\text{seed}}}$, we see that the shear flow dominated $x_{1/2}^{\prime }=0.145$. Plugging back in for $x^{\prime }$, we can solve for an expression for the seed value for ${\tau }_{\gamma }$ as follows:

These seed values are used as the initial guesses that are a required input of the lsqcurvefit function in MATLAB. An optional input of the lsqcurvefit function is a lower and upper bound for the fitting parameters which imposes bounds that these parameters cannot exceed. While one can imagine many experimental situations where biological or physical insight can provide useful bounds on these parameters, to remain unbiased and test the authenticity of our fitting models in a “worst case scenario,” initially we did not impose any bounds on the fitting parameters.

Simulation and fitting algorithms were run on a BlueHive supercomputer equipped with multiple CPU cores and nodes (Center for Integrated Research Computing, University of Rochester, Rochester, New York, United States). Utilizing high-throughput parallel computing allowed us to generate artificial MPFRAP data and fit to each of the three models for many different combinations of input diffusion coefficients, velocities, and shear rates. Simulations were run at 27 different flow speeds and 64 different shear rates, with 20 repetitions at each point, resulting in a total of 34,560 simulations. All software for simulating MPFRAP data and fitting that data are publicly available on Code Ocean: https://doi.org/10.24433/CO.4821156.v1.

2.4.

In vitro MPFRAP with Shear Flow

We also evaluated the three fitting models in an in vitro system with various combinations of velocity and shear rate. This system used a microfluidic channel with a height of $200\mu \mathrm{m}$, length of 50 mm, and width of 5 mm ($\mu$-slide I Luer, Ibidi). To initiate gravity-driven flow through the channel, a large inlet reservoir was placed on an adjustable lab jack to control the height of the reservoir relative to the channel. The outlet reservoir was kept at a constant height. The inlet reservoir was filled with 300 mL of a $1\mathrm{mg}/\mathrm{mL}$ solution of FITC conjugated to 2000 kDa dextran (Sigma). To this solution we added $1\mu \mathrm{L}/\mathrm{mL}$ red fluorescent microspheres (FluoSpheres, Molecular Probes/Invitrogen). The velocity of the bead/dye solution through the channel was manipulated by adjusting the height of the inlet reservoir and quantified by taking line scan images of the fluorescent beads, which were then analyzed in ImageJ. Three line scan images were taken prior to 3 MPFRAP measurements at a specific location within the channel, yielding at least 15 lines at each location. The measured velocities were plotted as a function of height within the channel and fit to a second-order polynomial, as seen in Fig. 3. The shear rate as a function of height within the channel was then calculated as the derivative of the flow profile. This produced a velocity and shear-rate at each location where a MPFRAP measurement was made, which was then converted to a scaled velocity ($v_s$) and scaled shear rate (${\gamma }_s$) using the focal volume dimensions of the lens.

Fig. 3

Experimental set-up for MPFRAP with shear flow and velocity and shear profile within a $200\mu \mathrm{m}$ microchannel. (a) Schematic diagram of experimental set-up for MPFRAP with shear flow. (b) Velocity measurements were computed from line scan images of fluorescent beads. Data are shown as mean ± standard deviation ($n\ge 15$). Velocity data were fit to a second-order polynomial and shear rate data were obtained by taking the absolute value of the derivative of the velocity profile.

A mode-locked Ti-Sapphire laser (Mai Tai; Spectra Physics, Mountain View, California, United States) delivering laser light with 80-fs pulses at a repetition rate of 100 MHz. Rapid modulation between bleach and monitor laser intensities was achieved with a KDP* Pockels Cell (model No. 350-80; Conoptics, Danbury, Connecticut, United States). Timing of the bleach and monitor pulses was delivered by a pulse generator (model No. DG535, Stanford Research Systems, Sunnyvale, California, United States). Simultaneously, the voltage output to the Pockels Cell was set and switched by a custom control box. After the Pockels Cell, the laser was directed through an Olympus Fluoview 300 laser-scanning microscope to the back aperture of a 0.8 NA, $40\times$ water immersion objective lens (Olympus, Center Valley, Pennsylvania, United States). The back aperture of the objective was overfilled, yielding a focal volume with $1/{\mathrm{e}}^2$ radius of $0.403\mu \mathrm{m}$ in the radial direction and $2.27\mu \mathrm{m}$ in the axial direction.¹⁰ The objective lens focused the focal volume at various positions within the microfluidic channel. The fluorescence emission was separated from the excitation light by a short-pass dichroic mirror (Chroma Technologies, Brattleboro, Vermont, United States). When imaging fluorescent beads to determine flow properties within the microfluidic channel, the emission signal was filtered by a 605/35 emission filter (Chroma Technologies, Brattleboro, Vermont, United States). When collecting FRAP data, the emission signal was filtered using a 534/30 emission filter (Chroma Technologies, Brattleboro, Vermont, United States) before being detected by a photomultiplier tube (PMT) (Hamamatsu, Bridgewater, New Jersey, United States). The PMT output was directed to a photon counter (model No. SR400; Stanford Research Systems, Sunnyvale, California, United States).

3.1.

Evaluating How the Diffusion-Convection Model Fails in the Presence of Shear

To evaluate the performance of the previous diffusion-convection model in the presence of shear, we performed a series of simulations modeling fluorescence recovery with various diffusion coefficients ranging from $0.5\mu {\mathrm{m}}^2/\mathrm{s}$ to $500\mu {\mathrm{m}}^2/\mathrm{s}$, with a flow velocity in the center of the focal volume of zero ($v_o=0$) and with a range of shear rates from 0.01 to $10000{\mathrm{s}}^{-1}$. These curves were fit to the diffusion-convection model and their diffusion coefficients were extracted. To determine accuracy of the fit, the fitted diffusion coefficient was divided by the input diffusion coefficient, thus an accurate fit will yield a ratio of 1. Figure 4(a) shows the accuracy of the diffusion-convection model for various diffusion coefficients as a function of shear rate. As expected, increasing shear will produce erroneously high values of measured D, to an extent that varies with the actual value of diffusion coefficient relative to shear rate. To remove the effect of diffusion coefficient, we scaled the $x$-axis by defining the scaled shear rate as the ratio of ${\tau }_D$ to ${\tau }_{\gamma }$, thus ${\gamma }_s=\gamma ({\omega }_r^2/8D)$. This dimensionless independent variable quantifies the relative contribution of shear and diffusion to the transport of molecules into and out of the focal volume to reveal a universal impact of shear rate on the accuracy of fit, which can be seen in Fig. 4(b). We see that the diffusion-convection model produces erroneous fits when ${\gamma }_s\sim 0.5$ and greater. Thus, when the recovery rate due to shear is more than half the recovery rate due to diffusion, the diffusion-convection model is no longer valid.

Fig. 4

Accuracy of diffusion-convection model for various diffusion coefficients as a function of shear rate and scaled shear rate. (a) MPFRAP curves were simulated with Monte Carlo simulation for various input diffusion coefficients with a central velocity of zero ($v_o=0$) for various ranges of shear rate. Curves were fit to the diffusion-convection model and the fit diffusion coefficient was extracted. Data are shown as mean ± standard deviation ($n=20$). (b) Scaling the $x$-axis of (a) by the relative contribution of shear and diffusion to the fluorescence recovery, ${\gamma }_s=\gamma ({\omega }_r^2/8D)$, allows curves of all values of diffusion coefficient to overlay into a single curve.

3.2.

Evaluating the Impact of Scaled Shear Rate on Previous Models

Using the concept of scaled shear rate, we moved on to observing how all three models behave as scaled shear rate is varied, at three different scaled velocities. The input diffusion coefficient for these simulations was set to $60\mu {\mathrm{m}}^2/\mathrm{s}$, which is a typical diffusion coefficient for FITC-BSA. Figure 4 shows the accuracy of fit for the diffusion-only model (red), diffusion-convection model (blue), and the new shear flow model (green) over a wide range of scaled shear values, for various scaled velocities, $v_s=0,0.5$, and 1. Similar to the scaled shear rate, the scaled velocity is a dimensionless independent variable that quantifies the relative contribution of flow and diffusion to the transport of molecules into and out of the focal volume¹⁰ (i.e., $v_s={\tau }_D/{\tau }_v={\omega }_r{v}_o/8D)$.

In Fig. 5(a), the scaled velocity is zero, representing an unphysical but mathematically attractive situation where the focal volume is placed exactly where the net transverse flow at $z=0$ is zero, reminiscent of a “whirlpool.” We see that all three models produce accurate diffusion coefficients up to a scaled shear rate of $\sim 0.3$. At scaled shear rates greater than $\sim 0.3$ and $\sim 0.5$, respectively, we see the diffusion-only and diffusion-convection model produce erroneous values of D, while the shear flow model continues producing accurate $D$ values up until ${\gamma }_s\sim 30$. In Fig. 5(b), the scaled velocity is set to 0.5 and we see the diffusion-only model produces erroneous values of $D$ at all values of scaled shear rate, while the diffusion-convection model produces accurate diffusion coefficients up until ${\gamma }_s\sim 0.5$, and the shear flow model produces accurate $D$ values up until ${\gamma }_s\sim 30$. In Fig. 5(c), the scaled velocity is set to 1, thus indicating an equal contribution of flow and diffusion to the fluorescence recovery in the center of the focal volume. Once again, we see the diffusion-only model producing erroneous $D$ at all values of scaled shear rate, while the diffusion-convection model produces accurate diffusion coefficients up until ${\gamma }_s\sim 0.5$, and the shear flow model continues to produce accurate diffusion coefficients up until ${\gamma }_s\sim 30$.

Fig. 5

Effect of scaled shear rate on accuracy of the three fitting models. Investigating the effects of scaled shear rate on the accuracy of fitted diffusion coefficient for the diffusion-only model (red, circles), diffusion-convection model (blue, triangles), and the shear flow model (green, squares) at $v_s=0$ (a), 0.5 (b), and 1.0 (c). Data are shown as the mean ± standard deviation ($n=20$).

3.3.

Accuracy of Fitting Models Over Large Ranges of $v_s$ and ${\gamma }_s$

Next, we performed Monte-Carlo simulations for a wide range of scaled velocities and scaled shear rates to visualize the area of the $v_s\text{-}{\gamma }_s$ parameter space over which the three models produced accurate diffusion coefficients. We preformed simulations at $27\times 64=1728$ different combinations of scaled velocity and scaled shear rate, respectively. We then fit the resultant fluorescence recovery curves to the three models and plotted the accuracy of the fitted diffusion coefficient as contour plots in Fig. 6. In these plots, using the $D_{\text{fit}}/D_{\text{in}}$ ratio as a measure of accuracy does not visually highlight erroneous values of $D_{\text{fit}}$ that are smaller than $D_{\text{in}}$ as well as it highlights $D_{\text{fit}}$ values that are too large. To facilitate visualization of the resultant accuracy of fit, here our measure of accuracy is presented as the square of the natural log of the fit diffusion coefficient divided by the input diffusion coefficient, $[\ln (D_{\text{fit}}/D_{\text{in}}){]}^2$. Thus, an accurate fit diffusion coefficient ($D_{\text{fit}}=D_{\text{in}}$) will yield a value of 0, if the fit diffusion coefficient is too small and the resulting $D_{\text{fit}}/D_{\text{in}}$ is 0.1, our metric will produce a value of 5.3, and if the fit diffusion coefficient is too large and the resulting $D_{\text{fit}}/D_{\text{in}}$ is 10, our metric will produce a value of 5.3. Thus, our metric displays erroneously small values of $D_{\text{fit}}$ in a similar way as erroneously large values of $D_{\text{fit}}$.

Fig. 6

Results of three fitting models over various scaled velocities and scaled shear rates. Fluorescence recovery curves were generated in a high-throughput manner at 1728 different combinations of scaled velocities and scaled shear rates ranging from ${10}^{-3}$ to ${10}^2$. Curves were then fit to the diffusion-only model (a, b), diffusion-convection model (c, d), and shear flow model (e, f), and the fitted diffusion coefficient was extracted. Data are shown as the mean of ($n=20$). Accuracy of the resultant diffusion coefficients is shown in panels (a, c, e) by the color bar, which represents ${[\ln (\frac{D_{\text{fit}}}{D_{\text{in}}})]}^2$, where an accurate fit diffusion coefficient will yield a value of 0, shown as dark blue. The average standard deviation is shown in panels (b, d, f) $(n=20)$, where smaller values are shown as dark blue.

Figure 6(a) shows the accuracy of $D_{\text{fit}}$ over the $v_s\text{-}{\gamma }_s$ parameter space when fit to the diffusion-only model. We can see that the diffusion-only model produces accurate diffusion coefficients for values of $v_s$ and ${\gamma }_s$ less than $\sim 0.3$ (dark blue region). Figure 6(c) show the accuracy of $D_{\text{fit}}$ over the $v_s\text{-}{\gamma }_s$ parameter space when fit to the diffusion-convection model. As expected, we see the range of scaled velocities that produce an accurate diffusion coefficient significantly increases over that of the diffusion-only model. However, for any value of $v_s$, the model produces erroneous diffusion coefficients when ${\gamma }_s$ exceeds $\sim 0.5$. Figure 6(e) shows the accuracy of $D_{\text{fit}}$ over the $v_s$-${\gamma }_s$ parameter space when fit to the shear flow model. We see that the shear flow model produces accurate $D_{\text{fit}}$ values over an area much larger than the diffusion-only and diffusion-convection models and begins to produce erroneous values when ${\gamma }_s$ exceeds $\sim 30$.

3.4.

Improving $D_{\text{fit}}$ with a priori Knowledge

As mentioned previously, no bounds were imposed on any fitting parameters in the above figures to remain unbiased and test the true authenticity and rigor of the shear flow model. Thus, Fig. 6(c) shows the accuracy of the shear flow model to determine the diffusion coefficient of a system in the worst case scenario where no a priori knowledge is provided. However, there is usually some level of a priori knowledge about one’s system that can be used to improve the fitting process. Table 1 shows how the fitted diffusion coefficient changes in the absence and presence of a priori knowledge about one’s system for two points taken from the red region of Fig. 6(e), where all three models failed to produce an accurate diffusion coefficient. First, we consider the case where we know within a factor of ten what the possible values are for all four fitting parameters. In this case, we impose a lower bound that is one tenth of the true values of each parameter and upper bound that is ten times the true values of each parameter. This results in a fitted diffusion coefficient that is much closer to the true diffusion coefficient ($D=60\mu {\mathrm{m}}^2/\mathrm{s}$). In the best case scenario, we suppose that the velocity and shear rate of a system is known or measured independently by the user (using fluorescent beads with line scans, etc.). In this case, ${\tau }_v$ and ${\tau }_{\gamma }$ can be fixed and ${\tau }_D$ and $\beta$ can be free fitting parameters with no bounds. We see that in this case we get a diffusion coefficient that is much more accurate than the diffusion coefficient obtained with no bounds on the fitting parameters.

Table 1

Fitted diffusion coefficient in the presence of different a priori knowledge about the system. Two regions of poor Dfit values were chosen from Fig. 6(c): one in with high vs and one with high γs. The first column shows Dfit when all four parameters are fit with no bounds [as in Fig. 5(c)], the second column shows Dfit when a lower bound (lb) and upper bound (ub) are imposed on all four parameters, and the third column shows Dfit when τv and τγ are known and τD and β are free fitting parameters with no bounds. The true diffusion coefficient is 60  μm2/s. Data are shown as mean ± SEM (n=20).

3.5.

In vitro MPFRAP with Shear Flow

Next, we performed MPFRAP in an in vitro system which allowed for various combinations of $v_s$ and ${\gamma }_s$ with a solution of 2000 kDa FITC-dextran with a known diffusion coefficient. The velocity of the solution through the channel was manipulated by adjusting the height of the inlet reservoir. To determine the shear rate at various locations within the channel, the velocity measured with line scan images was plotted as a function of height within the channel and fit to a second-order polynomial (Fig. 3); the shear rate as a function of height within the channel was computed as the derivative of the velocity profile, as shown in Fig. 3.

Multiple MPFRAP measurements were taken immediately after line scan acquisition at a specific height within the channel. The measured diffusion coefficient ($D_m$) was then compared to the known diffusion coefficient to determine the accuracy of the three MPFRAP models in vitro. First, we performed MPFRAP in a droplet of FITC-dextran to check that our system could accurately determine the diffusion coefficient in a diffusion-only system, this yielded a diffusion coefficient of $9.8\mu {\mathrm{m}}^2/\mathrm{s}$ which is comparable to values reported in the literature;²⁰ this measured diffusion coefficient ($D_{\text{truth}}$) was used as our ground truth value that all other measurements were compared to. All MPFRAP curves were fit to the three MPFRAP models with the same criteria as the simulated MPFRAP curves with no a priori bounds imposed, (representing the “worst case” scenario), as shown in Fig. 7.

Fig. 7

Accuracy of the three MPFRAP models in determining $D$ in vitro in the presence of various combinations of flow and shear rate. Scaled velocity and shear rates were measured as described in the methods section. Accuracy of $D_m$ is determined by calculating $\left[\ln \right(\frac{D_m}{D_{\text{truth}}}){]}^2$, which will yield values closer to zero when $D_m$ is closer to $D_{\text{truth}}$. Data are shown as mean ($n=3$) ± standard deviation. Note that bars for convective flow and shear flow in the first three columns are not visible because their mean values and error bars are so close to zero, while the bars for diffusion only and convective flow for the later columns are cut off, to keep the detailed behavior of the shear flow model visible.