Monte Carlo simulations of protein micropatterning in biomembranes: effects of immobile nanofeatures with reduced diffusivity

Nanoscopic features of reduced diffusivity have long been suggested to contribute to plasma membrane heterogeneity. Two prominent examples of this are highly dynamic lipid-mediated assemblies (‘membrane rafts’) and shells of annular lipids surrounding transmembrane proteins. Here, we simulated a micropatterning experiment, where such nanoscopic features are immobilized in specific areas within the live cell plasma membrane. We evaluated the effect of patterned nanofeatures of different sizes and diffusivities on the spatial distribution and two-dimensional mobility of tracer molecules. From this, we derive empirical models that describe the long-range tracer mobility as a function of the nanofeature density. In turn, our results facilitate the determination of nanofeature dimensions from micropatterning experiments.


Introduction
The perception of the plasma membrane has changed over the last decades from a fluid, rather homogeneous environment for lipids and proteins [1] towards a micro-and nanostructured matrix that controls cellular processes via the formation of spatial lipid/protein heterogeneities [2,3]. Various mechanisms were suggested to contribute to this compartmentalization, including the immobilization of membrane proteins via Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the cortical cytoskeleton [4] or exoskeleton fences [5], protein-protein interactions [6], lipid phase separation [7], or curvature-mediated formation of protein superstructures [8], to name a few.
Besides organizing the plasma membrane scaffold per se, non-homogeneous protein distribution further affects the mobility and spatial distribution of most other membrane constituents by defining diffusional barriers [9]; prominent examples include tight junctions in epithelia [10], and the axonal initial segment [11] as well as the synaptic membranes in neuronal cells [12]. In such protein meshes, the mobilityand hence permeability-becomes massively reduced for protein clusters compared to monomers [13], making such areas efficient size filters. Of note, transmembrane proteins are known to associate with an annular ring of lipids on a nanoscale [14,15]. Up to a hundred lipids may co-diffuse with the transmembrane protein forming a ring of reduced mobility extending several nanometers from the protein surface [16].
We have previously introduced live cell protein micropatterning as a tool to study the effective size of plasma membrane constituents-lipids as well as proteins-by determining the retardation of a fluorescently labelled tracer diffusing through regions of deliberately enriched and immobilized proteins [17,18]. In this approach, target proteins are rearranged in the live cell plasma membrane by growing cells on surfaces featuring micropatterns of antibodies. Micropatterned surfaces can be readily produced using soft lithography methods such as microcontact printing, reaching spatial resolutions down length scales of a few tens of nanometers [19], which allows for deliberate reorganization and enrichment of the target proteins to ON-regions, leaving OFF-regions depleted of target protein (figure 1(A)). Single particle tracking (SPT) of a fluorescently labelled tracer molecule can then be used to investigate tracer distribution and mobility in dependence of the nanofeature density.
In order to be able to extract target protein size as well as tracer-target interactions from diffusion data, we have previously employed extensive Monte Carlo simulations of tracers diffusing through regions of immobilized obstacles [20]. Here, we extended our previous study by addressing the effects of immobilized nanofeatures with reduced diffusivity, or lipid shells with reduced diffusivity co-immobilized with the target proteins ( figure 1(B)). We derived mobility ratios constants D ON /D OFF and tracer densities ρ ON /ρ OFF in ON-versus OFFregions and provide guidelines for the quantitative interpretation of such micropatterning experiments.

Simulation of nanofeatures
All simulations were performed in MATLAB (R2019a, The MathWorks Inc. Natick, MA) using a standard personal computer. Random walks of non-interacting point tracers were simulated on two-dimensional square maps containing random distributions of nanofeatures within the ON-region. The diffusion constant outside of the nanofeatures was set to D out (see table 1 for a list of variables). We simulated two different types of nanofeatures: (i) circles with diffusion constant D in < D out mimicking immobilized lipid rafts of radius R, (ii) circular rings (radius R) with diffusion constant D in < D out surrounding single central impermeable circular obstacles (radius R core ) mimicking immobilized transmembrane proteins including a ring of associated lipids ( figure 1(B)). In the following, index in indicates nanofeature segments of low diffusivity, index out membrane segments of high diffusivity (indicated in light grey and white in figure 1(B), respectively).

Implementation of the micropattern geometry
We considered here a micropatterning experiment, in which circular regions of immobilized nanofeatures are generated (termed ON-region) within regions devoid of such nanofeatures (OFF-region) ( figure 1(A)). In a real experiment, the size of the ON-region is on the order of a few micrometers, and hence more than hundred-fold larger than the immobilized nanofeatures. As output parameters of the simulation, we analyzed the ratio of tracer surface densities ρ ON /ρ OFF and the mobility ratio D ON /D OFF . For determination of ρ ON /ρ OFF , we simulated a single circular ON-region within a square OFFregion of identical area (dashed box in figure 1(A)). Note that in the absence of nanofeatures by definition the diffusion constant remains unaffected (D OFF = D out ); for determination of D ON /D OFF it was hence sufficient to simulate the diffusion behavior in ON-regions.

Distribution of nanofeatures
Nanofeatures were distributed with uniform probability in the ON-region. Overlapping of individual circular nanofeatures was permitted, yielding arbitrarily shaped patterns. In some cases, however, we were interested in the effects of clustered distributions of nanofeatures. To this end, nanofeatures with center coordinates (x, y) were positioned on randomly generated probability maps defined by p ( , where the number of clusters is denoted by N C , the uniformly distributed cluster centers by (x i , y i ), and the cluster size by σ. The density of nanofeatures was measured in terms of number density ρ or area fraction C. For practical reasons, we defined an intended area fraction C intend , which was used to calculate the simulated density of nanofeatures according to [21] Since this equation is only valid for uniformly distributed nanofeatures, we finally quantified the effectively covered area fraction C by creating pixelated masks with pixel size of 0.1 × 0.1 nm 2 , and counting pixels covered by the simulated nanofeatures. This was particularly relevant in case of clustered distributions of nanofeatures, where We calculated the total area fraction characterized by diffusion constant D in as C in . Additionally, we set C in in relation to the area fraction accessible to the tracer molecules 1 − C core ; we termed this ratio C ′ in . For simulation of type 1 nanofeatures mimicking lipid rafts, C ′ in = C in = C. In contrast, type 2 nanofeatures mimicking transmembrane proteins also contain an impermeable obstacle core with area fraction C core relates to the total area fraction C via C core = 1 − (1 − C) ) .

Simulation of tracer diffusion
We defined all lengths of our simulation in units of the nanofeature radius R. Diffusing tracer molecules were represented by point tracers with equal partitioning probability for membrane segments characterized by D in and D out . For all simulations we assumed periodic boundary conditions. At the beginning of each simulation, point tracers were distributed uniformly in the accessible space and started their random walk without mutual interactions. Off-lattice random walks were realized by repeatedly displacing tracers by a fixed step length l in a random direction, where l was chosen dependent on whether the tracer was located in a segment with D in , with D out , or crossing between the two. For each tracer position we determined whether the tracer was inside or outside of a nanofeature by calculating the distance d to the closest nanofeature center and comparing it to the nanofeature radius R. For d < R or d > R a tracer was considered inside or outside, respectively. Crossing was detected by a change in the tracer environment from inside to outside or vice versa. In order to adequately describe the transition between membrane segments of different diffusion constants, we assumed an underlying ballistic transport model. Briefly, the microscopic diffusion process can be interpreted as the result of a sub-nanoscopic ballistic motion in conjunction with elastic collisions on randomly distributed hypothetical particles. Hence, at this length scale l can be interpreted as the free path length of the tracer's ballistic motion. The density of these hypothetical particles defines l as well as the time for a single step ∆t. The step length and step time reflect the underlying sub-nanoscopic velocity v = l ∆t , yielding the microscopic diffusion constant D = l 2 4∆t = v 4 l = v 2 4 ∆t. Importantly, in this model the transition between segments of different diffusion constants is described without a change in velocity, but instead via a change in the density of the hypothetical particles. Hence, step lengths and times of two segments with different diffusion constants are related as Din Dout = lin lout = ∆tin ∆tout . Naturally, l out happens to represent the longest step length in the simulation, and was set to l out ⩽ R 10 . For tracers crossing from D out to D in , both the step length l cross and step time ∆t cross have to be modified. For this, we calculated the fraction f of a single step the tracer would have spent in the new environment, assuming no alteration of the tracer mobility at the border. Next, we resized step length and step time in order to account for the altered tracer mobility in the new environment, yielding l cross = (1 − f) · l out + f · l in and ∆t cross = (1 − f) · ∆t out + f · ∆t in (figure 2(A)). Crossing in the opposite direction was implemented analogously. Tracers detected inside a nanofeature core (i.e. d < R core ) were repositioned according to their ballistic reflection from the core surface.
In some cases, we were interested in the effects of spontaneously forming rafts. To this end, we employed a probabilistic trapping model [22]. Freely diffusing tracers were assumed to switch their mobility between D out and D in randomly in time, with characteristic transition time constants τ out and τ in , respectively. τ out and τ in were chosen to reflect average residence times as observed in simulations with immobilized nanofeatures. The time fraction β = τin τin+τout directly corresponds to the area fraction C covered by type 1 nanofeatures.

Analysis of simulated trajectories
In a real life SPT experiment, the described sub-nanoscopic motion of the tracers is not accessible. Instead, a typical experiment yields a sequence of positions of the tracer separated by time intervals t delay ≫ ∆t. We assumed here stroboscopic illumination with illumination times short enough to virtually freeze the motion of the particle during exposure. To the assumed diffusion model. lout and l in denote the mean free path length of the tracer's hypothetical sub-nanoscopic ballistic motion in the segments characterized by diffusion constants Dout and D in , respectively, ∆tout and ∆t in the according mean free time. The length of crossing steps at the interface between out-and in-segments was adapted as described in the main text. At the core surface (dark grey) tracers were assumed to be reflected. (B) The resulting trajectories were discretized, assuming stroboscopic illumination, using regular time intervals t delay , each containing at least 10 steps of length l in or lout (zoom in). For determination of macroscopic diffusion constants, we used t lag as described in the main text. mimic this experimental situation, we analyzed simulated trajectories by interrogating tracer positions ⇀ r (t) at regular time intervals t delay (figure 2(B)). Note that the differences between time steps ∆t in and ∆t out as well as the resizing of step times ∆t cross introduced small shifts in the discrete time basis for the different trajectories, preventing a precise definition of t delay . For convenience, we defined a set-point for t delay , selected the closest actual time-points within the trajectories, and took the average.
As common in SPT, square displacements were calculated for a large range of time lags with t lag = n · t delay . Mean square displacements were calculated as msd (t lag ) = r 2 (t lag ). Throughout this manuscript, we calculated the diffusion constant D ON by analysis of the mean square displacement as a function of the time lag t lag according to D ON = msdON 4t lag . Plotting the decadic logarithm of the mobility ratio D ON /D OFF as a function of log 10 (t lag ) resulted in sigmoidal curves (see Results ). The curves were approximated with the sigmoidal function dt, and a 1 , a 2 , a 3 , a 4 the free fit parameters. In this representation, physically interesting quantities can be derived from the fit parameters, such as the mobility ratio at low time lags, and the mobility ratio at large time lags We observed anomalous subdiffusion in the ON-regions In a double-logarithmic plot this proportionality translates to log 10 (D ON /D OFF ) ∝ (α − 1) · log 10 (t lag ). Consequently, the anomalous diffusion coefficient α can be directly extracted from the slope of this curve.

Used numbers
All simulations were carried out under periodic boundary conditions to reflect an infinite reservoir. The side length of the simulated square was set to 750 nm, which ensured that all observed results were insensitive to changes in the size of the simulated area. To emulate the size of large transmembrane proteins including a ring of annular lipids [16,23,24] or hypothesized lipid rafts [25], nanofeatures were simulated with a radius R = 2 nm to 15 nm. The tracer diffusion constant in absence of nanofeatures was set to D out = 0.4 µm 2 s −1 , a typical value for lipid diffusion in the plasma membrane of a living cell [26].

Results
In this manuscript, we used a Monte Carlo simulation approach to study the diffusional properties of membrane molecules moving in distinct membrane environments. This reflects a situation that can be experimentally realized in a micropatterning experiment. In such an experiment, a circular ON-region is defined by specific enrichment and immobilization of nanofeatures, whereas the surrounding OFF-region is left unchanged ( figure 1(A)). In a previous publication, we have analyzed the effect of immobilized obstacles on the diffusional properties of a tracer molecule focusing on steric hindrance and obstacle-tracer interactions [20]. Here, we went beyond this aspect and simulated nanofeatures, which are permeable to the tracer and slow down its diffusional motion. In the absence of any nanofeatures, the tracers were simulated with a diffusion constant D out . Nanofeatures were assumed either as circular segments with a reduced diffusion constant D in (type 1) or inert obstacles with a surrounding ring characterized by D in (type 2). If not noted otherwise, nanofeatures were distributed uniformly and allowed to overlap. Throughout this paper, we refer to the density of nanofeatures either in terms of number density ρ or in terms of area fraction C (see Methods).

Domains of reduced mobility
First, we were interested in the effects of immobilized type 1 nanofeatures on tracer mobility. A possible realization of this scenario would be immobilized lipid rafts in a micropatterning experiment. We simulated circular nanofeatures with a radius R = 8 nm and a reduced diffusion constant D in < D out . We performed simulations for increasing densities of nanofeatures and analyzed both the tracer surface density ratios ρ ON /ρ OFF and mobility ratios D ON /D OFF as a function of the nanofeature density C. We did not observe changes in the surface density ratio with increasing nanofeature density for any tested ratio of diffusion constants D in /D out (figure 3(A)). Since pure tracer retardation is insufficient for enrichment [27], this serves as validation of the ballistic transport model. Plotting the mobility ratio D ON /D OFF as a function of t lag , revealed anomalous subdiffusion according to D ON /D OFF ∝ t α−1 lag , with α < 1, over three orders of magnitude ( figure 3(B)). In the limits of very short and very long time lags, we observed apparent Brownian motion with α = 1. Fitting the curves for different nanofeature densities C with equation (3) yielded the plateau values of D ON /D OFF and the degree of anomalous diffusion α. We observed substantial deviation from Brownian motion with α down to 0.92 for large mobility differences between inside and outside of the nanofeatures (supplementary figures 1(A) and (B)).
In the limit of short lag times, tracers are limited to exploring their immediate surroundings typically without crossing nanofeature borders, yielding a plateau of pure Brownian motion with a mobility ratio D   (4) and (5)). Lines denote the empirical models for D ′ 0 (equation (6)) and D ′ ∞ (equation (7)).
constants D in and D out , characterized by the mobility ratio with D OFF = D out ( figure 3(C)). Given sufficient time, tracers can fully explore the heterogeneity of their environment, yielding apparent Brownian motion at reduced mobility ratio D ′ ∞ = lim t lag →∞ D ON /D OFF . Unfortunately, in a real life SPT experiment the complete time dependence of the diffusion constants is difficult to record, as very high spatiotemporal resolution of ∼1 nm and a few microseconds would be required to determine D ′ 0 . We were hence particularly interested in the dependence of the experimentally accessible parameter D ′ ∞ on the nanofeature density ( figure 3(D)). All simulated scenarios could be described as the sum of the weighted average D ′ 0 and an empirically derived parabolic correction term, yielding the mobility ratio The difference between D ′ ∞ and D ′ 0 can be interpreted as a measure of the spatial heterogeneity of the diffusion constant [28,29]. Varying the size of simulated membrane domains merely shifted the regions of anomalous subdiffusion in time without notable effect on D In a real-life experiment, micropatterned nanofeatures are often not homogeneously distributed across the areas of enrichment. To elucidate the effect of such imperfect patterns we simulated clustered distributions of nanofeatures ( figure  4(A)). Especially at low densities, the influence of clustering on the mobility ratio D ′ ∞ was negligible ( figure 4(B)). Note that for clustered nanofeatures at high densities ρ, the effective area fraction covered by nanofeatures differed significantly from C = 1 − exp ( −ρπR 2 ) due to increased overlapping of nanofeatures. Consequently, plotting the mobility ratio D ON /D OFF as a function of the number density ρ yielded  (7)). Arrow indicates the scenario shown in A. It turned out to be interesting to briefly explore the situation of spontaneously forming nanofeatures, which would correspond to dynamically forming lipid rafts. Again, we assumed these rafts to form exclusively in the ON-regions of the micropatterning experiment. To this end, we simulated heterogeneity in diffusion based on a probabilistic algorithm. Tracers were allowed to randomly switch between different mobilities D in and D out with characteristic transition times τ in and τ out , respectively. The fraction of tracers associated with a nanofeature is given by the time fraction β = τin τin+τout and directly corresponds to the area fraction C in the model of immobilized nanofeatures. Surprisingly, we did not observe any anomalous subdiffusion in this scenario, and the mobility ratio D ON /D OFF stayed at D ′ 0 for all simulated time lags (figure 5).

Inert obstacles with a surrounding ring of reduced mobility
On the nanoscale, transmembrane proteins associate with an annular shell of lipids, reducing the diffusivity around the protein up to several nanometers from the protein surface. To study the effects of immobilization and enrichment of such proteins in micropatterning experiments, we simulated ONregions with different densities of inert obstacles with radius R core featuring a surrounding ring of radius R with reduced diffusion constant D in < D out (type 2 nanofeatures) ( figure 1(B)). We tested several combinations of radii R core as well as mobility ratios D in /D out .
ρ ON /ρ OFF = 1 − C core (8) with denoting the area fraction covered by the solid obstacle cores (see Methods).
Plotting the mobility ratio D ON /D OFF as a function of t lag showed anomalous subdiffusion with α < 1. Particularly at high nanofeature densities, the degree of anomalous subdiffusion was higher than in the case of nanofeatures without a solid obstacle core (supplementary figures 1(C) and (D)). Considering the inaccessible area fraction C core , the mobility ratio at short time lags followed the weighted average with C ′ in = C−Ccore 1−Ccore denoting the fraction of accessible area covered by the viscous rings (see Methods) (figure 6(C)). With increasing density, tracer movements over long distances become increasingly impeded by impermeable obstacle cores until C core reached the percolation threshold C P = 0.676 [30] and long range paths were completely blocked (figure 6(D)). The mobility ratio at long time lags could be described by a combination of the models for inert obstacles and nanofeatures with reduced diffusivity according to the empirical equation .
(10) We also show the plateau values of the mobility ratio D ON /D OFF in dependence of the number density ρ to facilitate comparison with experimental results in supplementary figure 4.

Determining feature sizes with a model of inert obstacles
It turns out to be experimentally difficult to achieve nanofeature densities above C ≈ 0.4 [17,18]. Especially at low nanofeature densities, however, the ratio D in /D out and the nanofeature sizes R and R core have similar effects on the curves, making it difficult to discriminate different nanofeature types based on recording only D ′ ∞ . Without additional knowledge about the nanofeature morphology, a practical approach may be to plot D ′ ∞ as a function of ρ and fit the data with a model of inert circular obstacles, yielding the apparent obstacle size R app as a fit result. In the following, we attempted to analyze data for type 1 or type 2 nanofeatures by fitting with a model for inert obstacles.
In figure 7(A) we show the results for type 1 nanofeatures (here circles with radius R = 8 nm and diffusion constant D in = 0.1 · D out ; full line), compared to a model for inert obstacles (radius R = 7.94 nm; dashed line). Expectedly, considering typical experimental noise the two curves would be experimentally indistinguishable. In this case, the relative error Rapp−R R = −0.008 would be rather small. However, as shown in figure 7(B), the relative error may amount to much larger values for more faint differences in the diffusivity inside versus outside nanofeatures (black line). We also included type 2 nanofeatures (i.e. nanofeatures with an impermeable core): in this case, the inert obstacle model senses the effect of the core appropriately, yielding overall smaller deviations from the total size of the nanofeature.

Discussion
We provided here a first detailed theoretical evaluation of a hypothetical life cell experiment, in which viscous nanofeatures are immobilized at various surface densities in micrometer-sized patterns within the plasma membrane. We quantified the effect of different nanofeature morphology and size on tracer mobility and surface density, yielding the following key results: (i) As expected, the tracer density ratio ρ ON /ρ OFF is not affected by permeable nanofeatures (type 1), irrespective of the nanofeature diffusivity D in . This is in line with the common definition of a partition coefficient, which is independent of probe mobility [27]. In other words, in the absence of size exclusion or chemical attraction or repulsion, probe molecules would distribute homogenously over phase-separated membranes. (ii) Intuitively, one may have expected D ON /D OFF to be a weighted average of the diffusivities inside and outside the nanofeatures. While this holds true for the shortrange mobility ratio D ′ 0 , we found an interesting non-linear dependence of long-range mobility ratio D ′ ∞ on the nanofeature density (figures 3(C) and (D)). A similar phenomenon was previously reported in case of spatial heterogeneities in the diffusion constant [28,29]. The nonlinear dependence of D ′ ∞ on nanofeature density C is a consequence of anomalous subdiffusion in the presence of substantial heterogeneity in the diffusivity. This happens to be pronounced particularly at intermediate nanofeature densities (supplementary figure 1). Notably, this non-linear dependence vanished in case of spontaneously forming nanofeatures (figure 5), indicating that the stability of the spatial diffusivity map is a precondition for this phenomenon. (iii) We finally considered the situation of an experimentalist, who has recorded the ρ-dependence of D ′ ∞ without further knowledge about the underlying nanofeature morphology. According to Ockham's razor, the experimentalist may opt for fitting the simplest model, which would be the presence of immobilized impermeable obstacles, yielding the apparent obstacle radius as result. We found that under conditions of substantial diffusivity reduction inside the nanofeatures this value closely reflects the actual extension of the nanofeature (figure 7).
Our simulations are based on a few simplifying assumptions, which shall be discussed in the following.
(i) Consistently with our previous publication on the effects of inert immobile obstacles, we assumed point-like tracers throughout this manuscript. The experimental realization could be a single fluorescently labelled lipid molecule, which has a typical in-plane radius of less than 5 Å. (ii) While theoretically interesting, the assumption of immobilized type 1 nanofeatures may be experimentally difficult to achieve. Our underlying goal was to study the effects of raft immobilization on the mobility of lipid tracers. However, for experimental reasons the immobilization of the putative raft requires the presence of raft-resident markers, which can be captured on micropatterns. In most realizations, such a marker corresponds to a membrane protein, which constitutes an impermeable obstacle within the nanofeature; we termed this scenario nanofeature type 2. Still, in the case of a GPI-anchored protein this obstacle core can be fairly small (two fatty acid chains amounting to approximately 5 Å radius) compared to the nanofeature radius, hence closely resembling the situation of a type 1 nanofeature. (iii) We assumed here discrete transitions between D in and D out at the contour of the nanofeatures. This assumption simplifies putative fluctuations of the contours in case of lipid phase separation, or gradients in diffusivity in case of boundary lipids. While potentially having minor quantitative effects, we do not expect qualitative changes on the presented results. (iv) In this paper, we considered flat membranes. This appears justified at the bottom membrane of a cell anchored to the surface of a glass coverslip. Still, slight membrane undulations could be present, which would lead to a retardation of the apparent diffusion constant when measured via two-dimensional single particle tracking [31]. However, due to the high density of anchored nanofeatures in typical experimental settings, we do not expect substantial effects due to three-dimensional membrane topology. Ultimately, the problem could be approached by measuring the threedimensional diffusion, e.g. via supercritical angle fluorescence microscopy [32].
In conclusion, diffusion analysis in conjunction with a micropatterning experiment provides valuable information on the size of immobilized nanofeatures. Recording and including additional information such as tracer density ratios, and tracer diffusivities inside versus outside of the nanofeatures would further help to determine nanofeature dimensions more precisely.