Actin filaments partition primary cilia membranes into distinct fluid corrals

Cilia have been described as antennae that probe the physical environment surrounding cells and tissues (Pazour and Witman, 2003). To subserve this function, signal transduction proteins, including seven transmembrane-spanning G protein–coupled receptors (GPCRs) and their downstream effectors, are enriched within ciliary membranes, where they detect light, odorants, and secreted molecules, among other signals (Gerdes et al., 2009; Garcia-Gonzalo and Reiter, 2012; Pearring et al., 2013; Mukhopadhyay and Rohatgi, 2014; Reiter and Leroux, 2017). Signal transduction involves the diffusional encounter between membrane proteins. An activated GPCR sequentially collides with and activates G proteins, which then diffuse to modulate downstream effectors. This arrangement allows response amplification and provides sites for feedback regulation (Pugh and Lamb, 2000). The physical properties of cilia membranes are, therefore, central to understanding cilium signaling and function and cilia-related diseases, yet they remain poorly understood largely because primary cilia are smaller than the spatial resolution limit of conventional live-cell microscopy and because they are difficult to isolate to homogeneity.

Recent studies have taken advantage of the exquisite high spatial resolution afforded by single-molecule tracking to examine the transport of GPCRs and GPCR-like transmembrane proteins within cultured inner medullary collecting duct (IMCD) cell primary cilia and found that they mostly sample the ciliary membrane space via diffusion, with occasional, short-duration coupling to motor proteins (Ye et al., 2013; Milenkovic et al., 2015). It was also shown that ligand binding or receptor activation could result in local binding (Milenkovic et al., 2015) or coupling to intraflagellar transport (IFT) motors via BBsomes (Ye et al., 2018). More recently, a study of single IFT-88 molecules entering the cilium base and moving through the transition zone were found to undergo several changes in transport speed before emerging into the cilium, suggesting that several steps are involved in moving through the diffusion barrier at the cilium base (Yang et al., 2017 ,Preprint). The relatively slow tracking of proteins in these studies, however, did not allow analysis of the membrane structure or fluidity at a level previously achieved for the plasma membrane (Fujiwara et al., 2002; Kusumi et al., 2005, 2014; Suzuki et al., 2005).

Here we directly examined the physical properties of living IMCD3 cell primary cilia membranes using two independent, high-spatial-resolution and high-temporal-resolution approaches to examine the dynamics of intrinsic membrane proteins, tracking the movements of individual proteins within the ciliary membrane using quantum dots (Qdots) and a novel two-photon (2P) FRAP approach that we call 2P Super FRAP. We chose to examine the transport behavior of two-cilium–enriched GPCRs: the somatostatin receptor type 3 (SSTR3), a rhodopsin-like GPCR found within cilia of developing central nervous system neurons, where it plays a role in neural differentiation (Händel et al., 1999), and a chimera of the rhodopsin apoprotein, opsin (Rho), which is localized to photoreceptor sensory cilia (Palczewski, 2006; Deretic and Wang, 2012; Goldberg et al., 2016). Our results showed that intrinsic ciliary membrane proteins diffused rapidly within highly fluid local membrane domains delimited by actin filaments.

SSTR3-EGFP or a chimera of bovine rhodopsin with the third intracellular loop replaced by that of the SSTR3 receptor (Rhoi3S-EGFP; Geneva et al., 2017), each with an N-terminal myc tag (Fig. 1 A), were expressed in IMCD3 cells, and cilium elaboration was induced by serum starvation. The distribution of EGFP-labeled GPCRs appeared homogeneous in live ciliary membranes (Fig. 1 B), in agreement with previous studies (Geneva and Calvert, 2010; Hu et al., 2010; Trivedi et al., 2012; Ye et al., 2013; Geneva et al., 2017). GPCR dynamics were then probed using two independent approaches.

Several basic modes of GPCR transport along the ciliary membrane are possible, including Brownian motion, motor transport, and local binding to stationary objects. We began the study with an evaluation of the relative contributions of these modes to the overall movement of the GPCRs, tracking single molecules labeled with Qdots (Fig. 1, C–F). Tracking on our system allowed 20-nm precision (Fig. S1). Because the major modes of GPCR transport are most apparent in the axial direction, GPCR transport was examined in detail using axial kymographs (Fig. 1 D). SSTR3 and Rhoi3S displayed similar patterns of movement (Fig. S2, C and D; and Videos 3 and 4) that were distinct from motor protein–driven transport typified by components of the IFT system. The GPCRs interspersed Brownian movement with apparent processive or stalled movements of varying duration (Figs. 1 D and S2, C and D). IFT proteins generally traveled as focal puncta from one end of the cilium to the other without interruption (Fig. 1 G and Video 5). Examination of the instantaneous axial velocities of the GPCRs (Fig. 1 F) and IFT (Fig. 1 H) also showed significant differences. Anterograde and retrograde velocities of the GPCRs were symmetrical, peaking near zero and falling exponentially to a maximum rate of ~1–1.5 µm/s. In contrast, velocities of the motor-driven IFT proteins were approximately normally distributed, with mean velocities of 0.60 ± 0.02 µm/s (mean ± SEM) and 0.34 ± 0.02 µm/s, anterograde and retrograde, respectively, in good agreement with previous studies (Follit et al., 2006; Ye et al., 2013; Broekhuis et al., 2014).

MSD along the axial dimension of the cilium was performed to assess the major modes of GPCR transport. The analysis had qualitative predictions that gave some insight into the mechanisms of GPCR movement (Saxton and Jacobson, 1997; Kusumi et al., 2005, 2014; Metzler et al., 2014): If GPCR transport is mainly via diffusion, then the relationship between the time-averaged square displacement <x²> and the time step duration τ will be linear, with a slope proportional to the effective diffusion coefficient. If significant motor transport occurred, the relationship would initially arc upward, referred to as superdiffusion. Confinement of the GPCRs to local regions would result in a downward-arcing relation, known as subdiffusion. Invariably, the MSD(τ) plots for the individual GPCRs examined in this study had downward-arching shapes (Fig. 1, I and J). This result showed that the movements of the GPCRs along the ciliary membrane axis were “freer” at shorter τs and constrained at longer τs, and thus demonstrated subdiffusion-like behavior. Some degree of subdiffusion is expected, however, given the small size of primary cilia.

To examine the impact of boundaries at the base and tip of the cilia on GPCR motion, we calculated the expected MSD(τ) relation for diffusion within a cilium of finite length. The cilium membrane may be approximated as the surface of a right circular cylinder with diameter of ~300 nm. The chief variable in cilium membrane area, then, is the cilium length. The MSD(τ) relation for a membrane protein moving by diffusion along the ciliary axis x of length l_x may be described as

where σ = (2D_x)^0.5, D_x is the diffusion coefficient in the x direction, and

is the MSD(τ) at which the relation saturates at large τ because of the no-flux boundaries at the cilium base and tip, x = 0, and x = l_x (Kusumi et al., 1993). We compared the MSD(τ) predicted from Eq. 1 to the MSD(τ) obtained for each cilium (Fig. 1, I and J, green lines). The MSD(τ) relation of the tracking data followed the MSD(τ) predicted by Eq. 1 only through the first one or two τs (∼0.3–0.6 s), after which it peeled off to lower than predicted values, thus showing impediments to free axial diffusion other than the absolute length of the cilium (Table 1). Long-term GPCR tracking in most cases showed sampling of the entire cilium length (Fig. 1 K, fraction of cilium sampled [f_cs]; and Table 1). This indicates that the constraint on GPCR movements could not have resulted from absolute confinement to a subregion of the cilium.

A simulation of membrane protein transport on the cilium membrane was generated to examine the impact of various modes of GPCR transport on the MSD(τ) relation (Fig. 2 and Theory). The approach was based on a random walk on a cylinder surface into which other modes of transport, including transient local binding, transient motor transport, and transient confinement to corral domains, were incorporated. Results from the pure random walk simulation, based on the parameters of the SSTR3-expressing cilium displayed in Fig. 2 G, showed that the pure random walk contained many of the axial movement behaviors observed in tracking experiments (Fig. 2 B and Video 6), including apparently processive movements (Fig. 2 B, black arrows), some of which traversed the entire cilium length. Periods of low axial movements were also observed in the model output (Fig. 2 B, red arrows). Averaging ∼1,000 MSDs for a given τ produced results that largely agreed with the Eq. 1 prediction (Fig. 2 D), thus validating the random walk simulation. The results show that stretches of unidirectional movement and slow movement are inherent to the random walk and, alone, do not indicate discrete motor transport or local binding events.

Simulations that include transient GPCR association with motors were examined to understand the potential impact of active transport on the MSD(τ) plots. Motor displacement rate v_m = 0.6 µm/s was used, which represents the mean anterograde IFT velocity. Coupling of the GPCR to motors was allowed anywhere on the cilium surface with probability P_m = 0.2, the approximate frequency of motor coupling reported in Ye et al. (2013). The duration of the motor coupling, t_m, was varied from 1 to 12 time steps (0.3–3.6 s) and the model output was subjected to the MSD(τ) analysis. With the exception of t_m = 1, which produced an MSD(τ) plot that was indistinguishable from the pure Brownian motion prediction, all transient motor-driven simulations resulted in MSD(τ) that exceeded Eq. 1 prediction (Fig. 2 E and Video 7).

We explored the possibility that the subdiffusion observed in GPCR movement could be caused by transient local binding events. A model of stochastic, transient binding to immobile structures located anywhere on the cilium, with probability of binding P_b = 0.2 and binding duration t_b = 10 time steps (3 s), resulted in a nearly linear MSD(τ) relation (Fig. 2 F and Video 8). The apparent diffusion coefficient of the simulated GPCRs was uniform over all τ, a result that does not match the MSD(τ) behavior of the GPCRs in the living cilia.

To examine if there were infrequent motor coupling or local binding events not detected by MSD analysis, each tracking experiment was searched for clusters of sequential movements. We determined instantaneous velocities of the GPCRs and segregated them into two categories: those with velocities lower than the observed retrograde IFT rates (ν ≤ 0.2 µm/s) and those with higher velocities (Fig. 2 G). The durations and frequencies of sequential low-velocity events, and high-velocity events in either the anterograde or retrograde direction, were determined and histograms plotted (Fig. 2 H). The histograms of high-velocity clusters, in either direction, fell exponentially with cluster size, yielding a cluster decay constant, k_cl, of ∼1.2 time steps for both SSTR3 and Rhoi3S. Low-velocity cluster durations declined with k_cl = 1.3 time steps. Long-duration, low-velocity clusters were found, with 16 clusters of 10,000 having durations longer than 10 time steps and occasional low-velocity stretches of 36 time steps (∼11 s) or longer. This result indicates that GPCRs bind to relatively immobile structures with low frequency.

Remarkably, the temporal cluster analysis of the pure random walk model results did not differ substantially from GPCR tracking results (Fig. 2, I and J). The simulation results, however, did not contain the infrequent longer-duration, low-velocity clusters found in GPCR tracking, indicating that those events were likely genuine, in agreement with similar observations (Milenkovic et al., 2015). The simulations led to the conclusion that significant motor-based GPCR transport did not occur under our experimental conditions. This result stands for cilia expressing vastly different densities of the labeled GPCRs, from ∼100 to tens of thousands per cilium (Table 1), and thus does not appear to be an artifact of overexpression. The result is in agreement with a recent study that unliganded SSTR3 in IMCD3 cilia did not show coupling to motors (Ye et al., 2018).

Subdiffusion has been reported to arise from several mechanisms in live cells, including viscoelastic effects of membranes, local membrane curvature (Quemeneur et al., 2014), molecular crowding (Weiss et al., 2004), and corralling of the GPCRs into fluid membrane domains (Kusumi et al., 1993, 2014; Suzuki et al., 2005; Oh et al., 2014). With the exception of local corralling, the MSD(τ) relations for these mechanisms are expected to smoothly curve downward and to be well fitted with a single exponential with time exponent <1. Corralling into subdomains, on the other hand, is expected to result in MSD(τ) relations that exhibit more than one discrete phase, as originally shown using high-speed particle tracking (Kusumi et al., 1993; Oh et al., 2014).

We increased the sampling rate for Qdot tracking to 200 frames per second (fps), thus allowing τs down to 5 ms to be assessed (Fig. 3). At this frame rate, the Qdots appeared to stay within relatively small regions on the ciliary membranes, sampling them intensively, and then moving on to another discrete region (Fig. 3 B). The time spent within a given region varied from less than a second to about a second. Comparison of the apparent sizes of the regions to the positions recorded for the Qdot-labeled GPCR on a glutaraldehyde-fixed cell (Fig. 3 B, fixed) shows that the apparent confinement of the GPCRs was not caused by binding to stationary objects.

High-speed acquisition revealed two discrete phases of the MSD(τ) plot: a steep rise from 5 to ∼20 ms, followed by a dramatically shallower rise from ∼20 to 500 ms (Fig. 3 C). Linear regressions of the two phases showed that the first phases had, on average, 13-fold steeper slopes than the second phases, corresponding to mean diffusion coefficients of D_m1 = 2.9 ± 0.40 µm²/s and D_m2 = 0.24 ± 0.05 µm²/s (Fig. 3 E, n = 18; and Table 2). For this analysis, the MSD(τ) for Qdot-labeled GPCRs on glutaraldehyde-fixed cells was used to determine the noise floor of our imaging system and thus the lowest measurable MSD(τ) value (Fig. S1 D).

Previous studies of GPCR and lipid corralling on the plasma membrane estimated the corral size by direct observation of the tracking positions over time (Kusumi et al., 1993; Fujiwara et al., 2002; Suzuki et al., 2005). Taking this approach, examination of the clusters of GPCR positions in Fig. 3 B suggested that corrals may vary in size, with some clusters appearing quite compact and others more dispersed. However, several factors make drawing conclusions about corrals from direct observation of Qdot positions tenuous. The cylindrical geometry of the cilium makes it uncertain where the GPCRs are along the cilium circumference. The narrowness of the cilium reduced detectable corral jumping to one direction. Finally, and perhaps most importantly, as with the appearance of apparent processive motion, pure random walks may appear to have confinement regions, when observed by eye, when in fact they do not (Saxton, 1995).

We thus chose a statistical approach to detecting and estimating the mean corral size. We reasoned that, in the case of finite diffusion space, the MSD(τ) relation would reach a plateau at an MSD value that is proportional to the size of the diffusion space. This can be seen from Eq. 1, where the second term drops out at large τ:

and therefore

This result shows that the amplitude of the MSD(τ) plot at the break in slope, MSD_cor, is proportional to the axial length of the corral. On average, MSD_cor = 0.015 ± 0.002 µm², yielding an estimate of the mean size of the corrals of l_cor = 274 ± 20 nm (Table 2).

We next estimated the permeability of the corrals and the mean time that a GPCR remained within a given corral by fitting the high-temporal-resolution MSD(τ) plots with a confined random walk model (Fig. 4, Video 9, and Theory). We introduced a lattice of uniformly spaced corral boundaries along the length of the simulated ciliary membrane (Fig. 4 A), with period corresponding to the mean corral length, l_cor (Table 2). The ends of the cilium were treated as reflective no-flux boundaries, and the corral boundaries were treated as semipermeable, with the probability of GPCR crossing upon encounter being P_bc. P_bc and the diffusion coefficient for the GPCR within the corral, D_cor, were varied, and model MSD(τ) was compared with GPCR tracking results (Fig. 4 C). D_cor and P_bc were determined from root mean square error (RMSE) minimization (Fig. 4 D). The mean values from seven independent experiments were D_cor = 5.46 ± 0.44 µm²/s and P_bc = 2.8 × 10⁻³ ± 0.4 × 10⁻³ (Table 2). Importantly, the D_cor was slightly larger than D_m1 estimated by the initial slope of the high-temporal-resolution MSD(τ) analysis (Fig. 3 C and Table 2), indicating that even at 200 fps, the diffusion coefficient of the GPCRs within the corral was underestimated by the initial slope of the MSD(τ) plot.

Related to the probability of corral boundary crossing is the mean corral residence time, T_r. T_r was estimated from the simulations used to fit the high-temporal-resolution MSD(τ) plots by summing the time steps from corral entry, t_en, to corral exit, t_ex, and averaging over all corral visits. T_r varied broadly, as expected for stochastic boundary crossing (Fig. 4 B). On average, T_r = 0.71 ± 0.09 s (Table 2).

The Qdot tracking approach, which involves labeling of the GPCRs with a 20-nm-diameter Qdot coupled via two antibodies (Fig. 1 A), raises the possibility that extraciliary structures, including the glycocalyx, or receptor oligomerization may impact the results. To address these potential caveats, we used a completely independent, high-resolution approach to examining GPCR dynamics.

We noted that the mean corral width of ∼270 nm is on the order of the width of the 2P point spread function (psf) formed by focusing a titanium:sapphire (Ti:S) laser to the diffraction limit with a 1.2-NA water-immersion objective (Fig. 5, A and B; Calvert et al., 2007). Thus, exposure of EGFP-labeled GPCRs within the cilium to very brief, intense psf-shaped laser pulses would produce nonuniform, corral-sized bleaching patterns (Fig. 5 B). Thus, with sufficient sampling rates, ensemble GPCR movements within the corral domains could be directly monitored (see Materials and methods). The mean diffusion coefficient obtained with this approach, D Super FRAP or D_sf, was 1.86 ± 0.32 µm²/s (n = 7), in reasonable agreement with the fast phase from the single-molecule tracking results, D_m1 (Table 2).

We next examined the movement of ensembles of GPCRs on the ciliary membrane over longer distances and times. Here, we expressed SSTR3 labeled with the photoactivatable variant of GFP, PAGFP (Patterson and Lippincott-Schwartz, 2002) and the fluorescence recovery after photoactivation (FRAPa) approach we developed previously (Calvert et al., 2007). This high signal-to-noise FRAP approach allowed tracking of the GPCRs along the full length of the cilium, albeit at a much lower temporal resolution of 3 Hz. The time course of relaxation was fitted with the ensemble diffusion model with periodic semipermeable boundaries where the diffusion coefficients were reduced 13-fold (Fig. 5 I). The mean diffusion coefficient for SSTR3-PAGFP in this regimen was D_PA = 0.11 ± 0.013 µm²/s (n = 8), again in reasonable agreement with the estimate from phase 2 of the Qdot tracking results, D_m2 (Table 2). Thus, two independent approaches resulted in remarkably similar results at both short and long times and distances.

The results thus far are remarkably similar to those reported by Kusumi and colleagues (Kusumi et al., 2005; Ritchie et al., 2005; Suzuki et al., 2005; Hagiwara et al., 2011), who showed in a series of ground-breaking articles that lipids and the μ opioid receptor underwent hop diffusion between plasma membrane domains delimited by F-actin. Thus, we examined the possibility that the ciliary corrals were formed by actin filaments.

Treatment with 0.1 or 0.5 µM latrunculin A resulted in marked increases in the MSD(τ) slope, with MSD(τ = 4 s) significantly larger than in the pretreatment control in 7 of 11 trials (Fig. 6, A and G). Treatment with cytochalasin D gave a more nuanced result. 2 µM cytochalasin D on average had no effect (Fig. 6 G) where treatment with 13 µM resulted in significantly steeper MSD(τ) slope (Fig. 6, B and G). Treatment with 3 µM jasplakinolide invariably resulted in marked reduction in MSD(τ) slope and MSD(τ = 4 s) (Fig. 6, C and G). All treatments were compared with DMSO alone, which resulted in slight reduction in the MSD(τ) slope in two of six trials with no change in the remainder, for an overall mean of ∼10% reduction in mobility at MSD(τ = 4) (Fig. 6 G).

Fast tracking experiments showed that F-actin–perturbing drugs did not, however, change the corral size, as evident from the lack of change in the amplitude of the break between the two phases of the MSD(τ) plot (Fig. 6, D–F and H). Fitting of the MSD(τ) relations of the treated and nontreated cilia showed that the main impact of the drugs was to change the permeability of the corral boundaries; D_m1 was not changed (Fig. 6, I and J). Treatment of cilia with F-actin–perturbing drugs in FRAPa experiments resulted in significant changes in the effective diffusion coefficients (Fig. 7), thus independently confirming the tracking results.

These results strongly indicate that ciliary membrane corrals are delimited, at least in part, by F-actin and that the dynamics of F-actin may regulate the permeability of the corral boundaries.

It is not clear if F-actin is present in all cilia. To date, F-actin has been shown to be present in rod photoreceptor-connecting cilia (Chaitin et al., 1984) and within primary cilia under specific conditions where cilia shed ectosomes or their tips (Nager et al., 2017; Phua et al., 2017). F-actin–associated proteins were found to be present in primary cilia by a cilium membrane–specific proteomics approach, and actin itself was not (Kohli et al., 2017), although it was not shown if this membrane-based assay could detect actin associated with other membranes where it is known to reside. Phalloidin did not label cilia on fixed MDCK cells (Francis et al., 2011). However, to our knowledge, there are no studies that directly examine the presence of F-actin in living primary cilia under more physiological conditions. To determine if F-actin was present in the IMCD3 cilia in our study, we individually expressed two F-actin markers, Lifeact-mCherry and RFP-UtrCH, simultaneously with SSTR3-EGFP (Fig. 8 A). Both F-actin probes were detected in the live cilia by confocal microscopy (Fig. 8 A, white arrowheads).

However, because both of these probes are soluble when not bound to F-actin, and thus may simply enter the cilia by diffusion, as we have previously shown in photoreceptors (Najafi et al., 2012b), we confirmed that the ciliary signal indeed reflected F-actin binding. We compared the ratio of fluorescence in the cilia, F_cil, and in the cell body, F_cb, of the Lifeact-mCherry or RFP-UtrCH probes to that of unfused mKate2, which we used as a soluble protein control. Although weak mKate2 signal was also detected in the cilia (Fig. 8 A), the F_cil/F_cb ratio for Lifeact-mCherry or RFP-UtrCH was 20-fold higher than for that of mKate2 (Fig. 8 B), indicating specific F-actin binding.

We have shown that two intrinsic membrane proteins, the GPCRs SSTR3 and Rhoi3S, move within cilium membranes more quickly on short time and distance scales than over longer scales, which fits the broad definition of subdiffusion (Bouchaud and Georges, 1990; Feder et al., 1996; Saxton and Jacobson, 1997). A number of mechanisms may underlie subdiffusion of membrane proteins, including diffusion in an environment crowded by mobile or immobile obstacles (Saxton, 1987; Saxton and Jacobson, 1997; Banks and Fradin, 2005; Cho et al., 2012; Metzler et al., 2014) and viscoelastic effects of interactions with polymers, such as the cytoskeletal components (Weber et al., 2010). These mechanisms manifest MSD(τ) plots that are well described by a power law with respect to time step, τ^α, where α < 1 (Bouchaud and Georges, 1990; Feder et al., 1996; Saxton and Jacobson, 1997).

Our results with high-speed imaging show a strongly biphasic MSD(τ) that is clearly distinct from power law subdiffusion. This result is similar to previous studies that have reported corralled diffusion of membrane molecules on the plasma membrane. Fujiwara et al. (2002) showed that lipid molecules were transiently confined to 230-nm-diameter membrane regions that were themselves confined to 750-nm-diameter regions within the greater plasma membrane of rat kidney fibroblasts. Suzuki et al. (2005) demonstrated that the μ opioid receptor underwent “hop diffusion” between membrane corrals delimited by actin filaments. Later, Oh et al. (2014) showed that the serine chemoreceptor was transiently confined to a membrane region of 170-nm diameter that resided within a membrane compartment of 290-nm diameter within the polar regions of the Escherichia coli inner membrane.

Treatment of IMCD3 cells with the F-actin–perturbing drugs latrunculin A, cytochalasin D, and jasplakinolide showed that F-actin plays a major role in GPCR corralling in primary cilia. Moreover, the sizes of the corrals we detected, ∼270 nm, were similar to the corrals on the plasma membrane reported by Kusumi and colleagues (Kusumi et al., 2005; Ritchie et al., 2005; Suzuki et al., 2005; Hagiwara et al., 2011). Together with vastly different estimations for the effective diffusion coefficients obtained with high- and low-temporal-resolution 2P FRAP of multiple GPCRs, our results demonstrate that the primary cilium membrane possesses a lattice, much like that of the plasma membrane, that transiently corrals membrane proteins undergoing diffusion and that F-actin plays a significant mechanistic role in the corral structure. Corralled hop diffusion is thus a major mode of GPCR transport within IMCD3 cell primary cilia.

Whether or not F-actin is present within primary cilia is somewhat controversial. Attempts to stain for F-actin in ciliated cultured epithelial cells over the course of decades have failed to definitively show actin within cilia, or even within the periciliary membrane domain (Finck and Holtzer, 1961; Francis et al., 2011). However, a study of crane fly sperm flagella showed F-actin–like structures that bind heavy meromyosin (Forer and Behnke, 1972), and recent cilium proteomics studies showed that actin and actin-associated proteins are found within primary cilia (Konno et al., 2015; Kohli et al., 2017). Recent studies show that primary cilia appear to use actin filaments to shed ectosomes from their distal tips (Nager et al., 2017; Phua et al., 2017). Moreover, it is well established that photoreceptor-connecting cilia contain F-actin, where it appears to play a role in the elaboration of the disc membranes that contain the phototransduction cascade proteins (Chaitin et al., 1984; Chaitin and Bok, 1986; Williams et al., 1988; Chaitin, 1989; Nemet et al., 2014). Our results with mCherry-Lifeact and Tag-RFP-UtrCH showed that F-actin is present in the IMCD3 cilia and support its role as a regulator of corral boundary permeability.

A significant difference between plasma membrane and cilium membrane hop diffusion is the mean residence time, T_r, of the GPCRs within the corrals. In the plasma membrane, T_r for the μ opioid receptor was reported to be a mean of 45 ms for corrals ~210 nm wide and ~1.7 s for 730-nm-wide corrals, in which the smaller corrals were nested (Suzuki et al., 2005). Our estimate of T_r = 0.7 s in ciliary membranes is between these values. It is thus possible that further increasing the temporal resolution of our tracking system will resolve more than one corral structure. Additionally, although the actin destabilizing treatments we used resulted in increases in GPCR mobility at long times and distances, they did not result in increases in mobility that approached the high mobility estimated for the corral domain. The F-actin–destabilizing treatments were rather mild, and higher concentrations of the disrupting chemicals caused general cell morphology changes. Nevertheless, the fact that some mobility constraint remained after F-actin disruption suggests that additional factors may be involved in delimiting the ciliary membrane corrals. Supporting this idea, application of F-actin–perturbing drugs during high-speed tracking experiments showed that their impact was on the permeability of the corral boundaries; the apparent size of the corrals were unchanged (Fig. 6). This result suggests that F-actin may “gate” the corrals, rather than set the absolute coral size.

Interestingly, septin 2, the molecule proposed to form the diffusion barrier at the base of primary cilia (Hu et al., 2010), as well as septins 7 and 9, have recently been shown to be expressed along most of the length of the primary cilium (Ghossoub et al., 2013). Perhaps the septins contribute to corral boundaries. Transient GPCR confinement may also result from formation of ordered, nonfluid lipid raft domains of various sizes interspersed with less ordered, fluid membrane regions. Primary cilia are enriched in gangliosides GM₁ and GM₃ as well as the cholesterol binding protein prominin 1, suggesting that lipid rafts form in ciliary membranes (Janich and Corbeil, 2007). The observed nonsaturating subdiffusion could arise if GPCRs were excluded from the ordered domains and diffused only in the fluid regions around them, as was recently proposed in a theoretical paper on 2D diffusion in a poly-dispersed medium (Cho et al., 2012). Lipid rafts thus may contribute to corral formation or produce another layer of restricted diffusion beyond that of the F-actin described here. Further study is required to fully understand the ciliary corral boundaries at various time and distance scales.

Recent studies suggest that the primary cilium is divided into functional domains, perhaps to segregate signaling cascades, as evident from enrichment of subregions of the ciliary membrane with specific proteins. Inversin (Inv), a poly ankyrin repeat transcription factor involved in Wnt signaling, was shown to be localized to a region of the primary cilium distal to the transition zone and the ciliary necklace and proximal to the body and tip regions (Shiba et al., 2009). This localization pattern was the basis for defining an Inv domain in primary cilia where other proteins are also enriched, including Joubert syndrome–associated protein, ARL-13b/ARL-13, which were proposed to be localized to the Inv domain by a combination of active transport and leaky diffusion barriers (Cevik et al., 2013). These and similar results have led to the proposal that the primary cilium contains at least five distinct subdomains with different functions and different complements of enriched proteins (Blacque and Sanders, 2014). How these domains are established is not currently known; it would be of interest to see if F-actin–perturbing drugs disrupt these domains.

One of the remarkable features of the ciliary corrals is the fluidity of the corralled membrane domains. The mean diffusion coefficient of the GPCRs at short times and distances obtained with high-speed tracking, D_cor = 5.5 µm²/s, is among the largest diffusion coefficients reported for intrinsic membrane proteins. It is more than 10-fold larger than the diffusion coefficient for rhodopsin within photoreceptor disc membranes (Liebman and Entine, 1974; Poo and Cone, 1974; Wey et al., 1981; Gupta and Williams, 1990; Wang et al., 2008; Govardovskii et al., 2009; Najafi et al., 2012a). Similar membrane diffusion coefficients, averaging ∼4.2 µm²/s and ranging up to 30 µm²/s or more, were reported for the μ opioid receptor during high-speed single-molecule tracking in the plasma membranes of normal rat kidney fibroblast and Chinese hamster ovary cells (Suzuki et al., 2005). Given the large diffusion coefficient for the GPCRs in our study, the mean residence time of a GPCR within a corral and the mean corral size, we predicted that the corral membrane area fraction sampled by the GPCRs during a corral visit was ∼0.7 (Fig. S3). This shows that confinement of the GPCRs within corrals will dramatically increase the likelihood of a given GPCR to encounter other proteins, such as a G protein, within the same corral, and would reduce the probability of encountering one outside of the corral. Thus, the corrals may help to generate local signaling nanodomains and could conceivably play regulatory roles in ciliary signaling by limiting the number of interactions between cascade components on the time scale of signal transduction.

Protein expression of all constructs was driven by the CMV promoter. The SSTR3-EGFP fusion construct was a gift from K. Mykytyn (Ohio State University, Columbus, OH). The SSTR3-mKate2 fusion protein was created by subcloning mKate2, a gift from E.L. Snapp (Albert Einstein College of Medicine, New York, NY), in place of EGFP in the SSTR3-EGFP construct. The bovine rhodopsin DNA construct with EGFP fused to rhodopsin with the last eight amino acids of rhodopsin appended to the C terminus of EGFP (Trivedi and Williams, 2010) was obtained from D.S. Williams (University of California, Los Angeles, School of Medicine, Los Angeles, CA). A rhodopsin chimera with third intracellular loop of SSTR3 (T236-R258) replacing the third intracellular loop of rhodopsin (E232-K248), Rho-i3S, was constructed using the FastCloning method (Li et al., 2011). Replacing the i3 loop on rhodopsin with the i3 loop pf SSTR3 was necessary to achieve Rho-EGFP enrichment in the IMCD3 cilia. For single-molecule tracking, the SSTR3-EGFP and Rho-i3S chimera were subcloned into PrecisionShuttle mammalian vector with N-terminal Myc tag (PS100012; OriGene). To mark the base of the primary cilia, GFP-xcent (a gift from J. Amack, SUNY Upstate Medical University, Syracuse, NY) and pTagBFP-C (a gift from A. Gross, University of Alabama, Birmingham, AL) was subcloned to make BFP-xcent construct. IFT20-GFP and IFT88-GFP constructs were gifts from G. Pazour (University of Massachusetts Medical School, Worcester, MA). pCS2.0 Tag-RFP-T-UtrCH was obtained from Addgene (plasmid 101279). mCherry-Lifeact-7 was obtained from Addgene (plasmid 54491).

IMCD3 cells were grown in DMEM/F12 culture medium with 10% FBS, penicillin/streptomycin, l-glutamine, sodium pyruvate, and nonessential amino acid solution at 37°C. Transfection was achieved with Lipofectamine 2000 (Life Technologies) or jetPRIME (Polypus Transfection). Before imaging, cells were seeded onto Fisherbrand No. 1, 25-mm cover glass (12-545-86), grown to confluence, and then starved by reducing FBS to 0.2% for 12–48 h to induce cilium elaboration.

Cells were incubated for 30 min in blocking buffer (DMEM/F12 1:1, with l-glutamine, Hepes, without Phenol Red, 1% BSA); incubated for 15 min in blocking buffer with 1:20,000 dilution of mouse anti–c-myc antibody (Covance); rinsed three times with blocking buffer; incubated for 10 min with 1:12,000 Qdot 625-nm goat F(ab′)2 anti-mouse IgG (A10195; Thermo Fisher); and washed three times. This protocol achieved sparse GPCR labeling. The bases of cilia were identified just before Qdot tracking by taking images of the cilium in the BFP (BFP-xcent) and EGFP (GPCR) channels. Cilia that lay flat relative to the imaging plane were selected for analysis. Tracking of Qdot 625 was performed at 37°C with an inverted epifluorescence microscope (Olympus IX70) using a 60× 1.45-NA Plan Apo objective. The ∼1-µm depth of focus means that cilia were within ∼10 degrees of the z-plane. Excitation was broadband in a window between 400 and 490 nm. The frame rate was set to 3.333 fps (image acquired every 0.3 s), 5 fps (image acquired every 0.2 s), or 200 fps (image acquired every 5 ms). The low frame rates allowed use of lower excitation intensities: excitation power was <3 mW, because of longer image integration times, thus minimizing photobleaching of the EGFP label that identified cilium position. This allowed tracking over tens of minutes and was an appropriate frame rate for assessing motor-driven movements or local binding events. High-frame-rate (200-fps) acquisition required at least 10-fold higher excitation intensities to offset the lower image integration times. This resulted in more rapid EGFP bleaching and thus prevented tracking for longer than ∼30 s. Images of the red (Qdot) and green (EGFP) channels were acquired simultaneously with an Optosplit III (89North) with a custom set of dichroic mirrors and barrier filters that were arranged to simultaneously project the red and green channels onto different regions of an EMCCD camera chip (iXon DU-897; Andor) running on Nikon Elements software.

After tracking acquisition, the images were converted to tiff format, and the fluorescent wavelength channels that were recorded simultaneously on the EMCCD chip were separated into individual files using the Fiji Cairn Image Splitter plugin. The Qdot positions in each time series image were determined and linked together to form contiguous trajectories using a particle tracking code written in Matlab (Crocker and Grier, 1996). The Qdot positions were then layered onto the images from the EGFP channel, which shows the position of the cilium, and the displacements along the ciliary axis were determined using a custom Matlab code developed in our laboratory.

The MSD as a function of time steps τ was calculated according to

where τ is the time step interval for the position difference, x is the axial position of the Qdot or the simulated GPCR (in the case of the model), N = n − τ is the number of positional differences, and n is the total number of image time steps. This approach averages all pairs of points with a given time step interval τ (Lee et al., 1991; Zhang et al., 1991; Saxton and Jacobson, 1997).

Live confocal imaging of mCherry-Lifeact and Tag-RFP-UtrCH was performed on an LSM510 using a Plan-Apochromat 63×/1.4-NA oil-immersion objective (Zeiss) at 30°C. All other confocal and 2P imaging was performed with a custom-built confocal/multiphoton microscope (Peet et al., 2004; Calvert et al., 2007, 2010) consisting of an argon-ion laser (model 163C; Newport Corp.) and a Ti:S laser (Mai Tai HP; Newport Corp.), each of which were focused to the diffraction limit with a 60×, 1.2-NA, water-immersion objective (Plan Apo VC; Nikon). Fluorescence emission was detected with photon-counting avalanche photodiodes (SPCM-AQR-14; PerkinElmer). Data acquisition and instrument control were implemented with computer acquisition and control boards and custom software written in LabView (National Instruments) in collaboration with M. Coleman (Coleman Technologies, Newtown Square, PA). Cells were maintained at 37°C. For FRAP analyses, the 3D coordinates for the photoconversion pulse were manually selected from initial 3D scans of the cells using the 488-nm line of the argon ion laser. 2P laser exposures were controlled by a high-voltage, fast Pockels cell (350-50; Conoptics) capable of modulating the beam power on the microsecond time scale.

The 2P Super FRAP approach is a modification of one we previously described, the point-blast-point-scan FRAP method (Calvert et al., 2007), and is similar to early photobleaching approaches to membrane protein dynamics using rhodopsin optical density (Poo and Cone, 1973, 1974; Liebman and Entine, 1974) or fluorescence photobleaching recovery (Axelrod et al., 1976a,b), with the exception that our approach has subdiffraction-limit spatial resolution and very high temporal resolution, achieving 200-kHz acquisition rate.

Cilia expressing SSTR3-EGFP were scanned in 3D in confocal mode, and a point along the cilium length was selected in x,y,z for 2P Super FRAP (Fig. 5 C). Fluorescence was monitored at that point using the Ti:S laser tuned to 920 nm, ∼3 mW, with an acquisition rate of 200 kHz. A series of 25–50 bleach pulses of 0.1–10-ms duration, ∼30 mW, were applied at time 0, and the recorded traces were normalized to the prebleach fluorescence and averaged (Fig. 5 D, black trace). To reduce oversampling and noise, the data were decimated: i.e., smoothened by 10-window boxcar averaging and then 10-fold downsampled (Fig. 5, D and E, symbols). Minimal diffusion of SSTR3-EGFP occurred during the short bleach exposures, thus allowing fluorescence changes to be monitored as the psf intensity profile bleach pattern relaxed into an approximately homogeneous distribution. The recovery time course was then fitted with the numerical ensemble diffusion model (Fig. 5, E and H; and Theory).

FRAPa was implemented using the point-blast region of interest scan method (Najafi et al., 2012a). PAGFP-labeled SSTR3 was found in cilia by 3D confocal scanning using the 488-nm line of an argon-ion laser. A position along the cilium was selected in 3D where SSTR3-PAGFP was photoconverted from the minimally fluorescent state to the fluorescent state by 0.1–10-ms, 10-mW mean power pulses from the Ti:S laser tuned to 820 nm. The equilibration of the activated molecules was then monitored with serial x–y confocal scans using 488-nm excitation. Rapid focus corrections were made before and after the photoconversion pulse to account for the measured focus difference between 820- and 488-nm illumination caused by objective chromatic aberration (Calvert et al., 2007).

Raw images were processed using custom Matlab (MathWorks) routines to correct for slight field inhomogeneities inherent to the optical system and for nonlinearities in the photon detectors as follows. Full-field scans of solutions of fluorescein were averaged to generate image field-flattening maps. Spatial photon counts in the maps were normalized to the maximum counts, and acquired images were multiplied by the maps (Calvert et al., 2007, 2010; Najafi et al., 2012a). No other image processing was done.

Ensemble diffusion model solutions were computed for various values of D on an appropriate time grid and convolved in 3D with the measured multiphoton psf of the Ti:S laser (2P Super FRAP experiments) or the measured psf of the argon-ion (488-nm) laser (PAGFP FRAPa experiments) to predict the measured fluorescence relaxation. The RMSEs between data and theory were calculated as follows: nine data samples centered on the time at which each theory point was computed were averaged, and the difference between theory and the averaged data point was computed, squared, summed, and divided by the total number of theory points. The optimal fitted trace was found as that which minimized the RMSE (Calvert et al., 2007).

For 2P Super FRAP and PAGFP FRAPa experiments, only SSTR3-EGFP– or SSTR3-PAGFP–expressing cells were analyzed because they more efficiently enrich in the cilia, leaving virtually no signal in the apical membrane and other cellular structures (Geneva et al., 2017).

Direct measurement of the psf profile in our confocal microscope was performed by scanning 0.1-µm-diameter fluorescent microspheres. Fluorescence profiles of the microspheres were fitted with Gaussians, yielding standard deviations (σ) in xy of σ_xy = 0.14 µm and in z of σ_z = 0.56 µm. To estimate the EGFP or PAGFP photoconversion profile on the cilium membrane (Fig. 5), the psf intensity profile was approximated as a 3D Gaussian:

where a is the intensity maximum at the 3D center of the profile. For this analysis, we used a normalized psf profile where

The estimate of the profile on the cilium was then obtained from the interpolation of the psf at the surface of a right circular cylinder with diameter of 300 nm whose axis was perpendicular to the z-axis of the psf.

Results are presented as mean ± SEM. When pairwise comparisons were made, significance was tested with appropriate t tests. When multiple means were compared, a one-way analysis of variance was performed with α = 0.05, where the null hypothesis was that all mean values were equal. Upon rejection of the null hypothesis with P < 0.05, post hoc comparisons among the mean values were performed by two-sample t test. To account for the variance among the multiple pairs of mean values compared, the Bonferroni correction was applied to α such that α_b = α/m, where α_b is the corrected α and m is the number of comparisons being made among the means.

We used two high-spatiotemporal-resolution approaches to explore the dynamics of intrinsic membrane proteins on live primary cilia, single-particle tracking and 2P Super FRAP. In this section, we describe the models of protein transport we used to quantitatively evaluate the results and differentiate between potential modes of transport.

The primary cilium membrane is essentially the surface of a right circular cylinder and thus is conveniently described in cylindrical coordinates. The spatial coordinate system is defined as $r ⇀ =( r,z,θ ),$ where r is the radial distance, z is the axial distance, and θ is the angle about the axis. Owing to the fact that GPCRs are intrinsic membrane proteins and thus confined to the ciliary membrane at r = r_m, the radius of the ciliary membrane, modeling of GPCR transport is reduced to two spatial dimensions where z and θ are the only spatial variables.

The simulation was based on a 2D, continuous random walk on the surface of a right circular cylinder and implemented in Matlab or the C programming languages (Fig. 2 A). A point on the cylindrical surface, p(r_m,θ,z), was randomly chosen as the initial position of a simulated GPCR, and the direction of movement away from the initial and all subsequent positions was determined by a randomly selected angle, φ. The simulated GPCR was then moved a distance d along the φ trajectory on the cylinder surface as determined by d = vt_s, where t_s is the magnitude of the time step and v is the velocity of the GPCR. v was defined by the estimate of the membrane diffusion coefficient D_m1 as:

The cylinder ends, corresponding to the base and tip of the cilium, were treated as reflective, no-flux boundaries. Axial positions of the simulated GPCR (Fig. 2 B) were thus equivalent to the axial kymographs of the GPCR tracking shown in Fig. 1.

To examine the impact of transient motor-driven transport, confinement or local binding on the predicted MSD(τ) relation, stochastic coupling to motors, binding to immobile objects, or transient confinement to axial subdomains were incorporated into the random walk simulation. To accomplish the addition of a motor-driven component, coupling to an anterograde or retrograde motor occurred stochastically with a probability P_m applied at each time step that was arrived at by free diffusion. Once on a motor, the simulated GPCR moved axially a distance d_m = ±v_m × t_m, where v_m is the motor displacement velocity. The simulated GPCR moved in a + or – direction for the number of sequential time steps defined by t_m, the motor coupling time.

The impact of transient local binding was examined by implementing stochastic coupling of the GPCRs to immobile objects. Here, a freely diffusing molecule executing the random walk bound to an immobile object at arbitrary position within the cilium membrane with a probability P_b. Once bound, the molecule remained stationary for a duration t_b after which it resumed the random walk until it again became bound.

We examined the impact of transient corralled diffusion by generating a theoretical ciliary membrane onto which we imposed corral boundaries at arbitrary spacing and with specified permeability to the GPCRs. The boundaries were approximately square, where sides parallel to the axis of the cylinder had a length of L₁, and sides perpendicular to the axis had a length of L₂ = C/Round(C/L₁). Simulated GPCRs executed random walks within the corral domain as described above. When the path of the GPCR movement crossed a boundary, a random decision was made as to whether the GPCR was reflected back into the same corral in which it was diffusing or crossed the boundary and entered the neighboring corral. The boundary reflection versus crossing decision was governed by a probability of corral boundary crossing, P_bc. Matlab code for the random walk simulation is provided in the online supplemental material.

The model of ensemble diffusion on the cilium surface was based on molecular diffusion models we previously published (Calvert et al., 2007, 2010; Najafi et al., 2012a), where we described diffusion in 3D in a cylinder. In the present work, diffusion is limited to two spatial dimensions at the cilium surface, and reduced flux corral boundaries were introduced. The general diffusion equation may be written as

where c is the concentration of diffusing substance, t is time, D is the diffusion coefficient, and $∇2$ is the Laplace operator. Eq. 2 is a parabolic partial differential equation that may be applied in any coordinate system. The general diffusion equation in cylindrical coordinates may be written

where D_r, D_θ, and D_z are the diffusion coefficients in the indicated directions, and Q_s is a source or sink of the diffusing substance. Because the GPCRs are intrinsic membrane proteins, and thus are confined to the ciliary membrane, D_r = 0 and the GPCRs diffuse only in z and θ at r = r_m, the radius of the cylindrical cilium at the membrane. Eq. 3 thus reduces to

Eq. 4 is first order in time and second order in space, requiring one initial condition (IC) and two boundary conditions (BC) in each spatial variable. The IC is

The BCs in θ are

which reflect the symmetry in the source term and assume isotropic diffusion in θ and z. The BCs in z are

which specify no-flux, Neumann boundary conditions at the base and tip of the cilium. The assumption of no-flux at the cilium base appears appropriate on the time frame of the experiments described in light of the presence of a significant diffusional impediment for membrane proteins between the ciliary membrane and the larger apical membrane (Geneva and Calvert, 2010; Hu et al., 2010).

Solutions of Eqs. 4–7 were found using the numerical method of lines, where the partial differential equations specified in Eq. 4 were replaced with ordinary differential equations (ODEs), finite differences approximations of the spatial derivatives (Schiesser, 1991; Schiesser and Griffiths, 2009). For the general grid point p(i,j), where i and j are indices for θ and z, respectively, the spatial derivative approximations were

and

where Δθ and Δz are the spatial increments in θ and z. Note that in the present treatment, the inherent singularities at r = 0 were avoided because of the absence of radial diffusion and because r_m ≠ r₀. For the complete numerical solution, which includes radial diffusion and the regularization of those singularities, see Calvert et al. (2010).

Approximations of the grid points neighboring the no-flux boundaries at the base and tip of the cilium were

at the cylinder base, where the flux between (i,j – 1), a point that lies outside of the base of the cilium, and (i,j) is taken to be 0, and

at the cylinder tip, where the flux between (i,j + 1), a point that lies outside of the tip of the cilium, and (i,j) is taken to be 0. The inhomogeneous source term was

where $α m *$ is the effective 2P cross section of the fluorophore, γ* is the quantum efficiency of photoconversion, and I(i,j,r = r_m) is the photon flux density at the cilium membrane. Substituting Eqs. 8–12 into Eq. 4 gives a system of n_θ × n_z ODEs in t. The ODEs were programmed in Matlab and integrated forward in t from t = 0 using the Matlab ODE integrators ode15s or ode45.

The model described above is appropriate for diffusion on the ciliary membrane in the absence of corrals. To assess the impact of corrals on ensemble diffusion, we mapped boundaries onto the cylindrical grid that acted to reduce the local flux of theoretical GPCRs. We mapped approximately square corral boundaries onto the cylinder surface where sides parallel to the axis of the cylinder had a length of L₁, and sides perpendicular to the axis had a length of L₂ = C/Round(C/L₁). Transport across the boundaries was diffusive, with diffusion coefficient D_cb distinct from the bulk membrane diffusion constant D_m. Solutions to the diffusion equations with local change in diffusivity at the boundaries were found numerically, as follows. The flux at the grid point neighboring a corral boundary in the positive z direction was approximated as

and at the grid point neighboring the corral boundary in the negative z direction was approximated as

The flux at the grid point neighboring the corral boundary in the positive θ direction was approximated as

and the flux at the grid point neighboring the corral boundary in the negative θ direction was approximated as

Flux between grid points within the corral boundaries was found by solving Eqs. 4, 13, and 14 as described above, where D_z and D_θ were replaced with D_cb. Matlab code for the enesmble diffusion model is provided in the online supplemental material.

Fig. S1 shows positional precision of the tracking system. In Fig. S2, axial kymographs of GPCR tracking results reveal a mixture of Brownian motion, apparent processivity, and apparent binding to immobile structures. Fig. S3 shows impact of corrals on sampling of the ciliary membranes by GPCRs. Video 1 is a high-magnification video of the movement of Rhoi3S labeled with a Qdot625. Video 2 is the same as Video 1 with tracking overlaid. Video 3 is a low-magnification video of SSTR3 labeled with Qdot625. Video 4 is a low-magnification video of Rhoi3S labeled with Qdot625. Video 5 shows tracking of IFT20-GFP. Video 6 shows output of a random walk simulation. Video 7 shows output of a random walk simulation with stochastic motor coupling. Video 8 shows output of a random walk simulation with stochastic binding to immobile objects. Video 9 shows output of a random walk simulation with semipermeable corrals. Matlab code for the random walk and ensemble diffusion models are also provided as a ZIP file.

We thank members of the Center for Vision Research and Mira Krendel for helpful discussions and Zach Fechtner for technical assistance.

This work was supported by National Institutes of Health grant R01EY018421 (P.D. Calvert), a Postdoctoral Award from Fight for Sight (H.Y. Tan), and an unrestricted grant to the State University of New York Upstate Department of Ophthalmology from Research to Prevent Blindness. We are grateful for support from Lions Clubs of Central New York, District 20-Y1.

The authors declare no competing financial interests.

Author contributions: Conceptualization: P.D. Calvert; methodology: S. Lee, H.Y. Tan, P.D. Calvert; software: S. Lee, H.Y. Tan, A. Kruglov, P.D. Calvert; validation: S. Lee, H.Y. Tan, P.D. Calvert; formal analysis: S. Lee, H.Y. Tan, I.I. Geneva, A. Kruglov, P.D. Calvert; investigation: S. Lee, H.Y. Tan, I.I. Geneva, P.D. Calvert; resources: P.D. Calvert; data curation: P.D. Calvert; writing-original draft: S. Lee, P.D. Calvert; writing – review and editing, S. Lee, P.D. Calvert; visualization: S. Lee, P.D. Calvert; supervision: P.D. Calvert; project administration: P.D. Calvert; funding acquisition: P.D. Calvert.