Interacting passive advective scalars in an active medium
Sk Raj Hossein Rituparno Mandal Madan Rao
physics.bio-ph cond-mat.soft












    n J∂/ł()̊∇⃗𝐫̂𝐞̂_θ
            ∂/⟩⟨αβ̱∇̣∂ nττ_vτ_aτ_dv̌ x f Q D J V∑_α∑_αβ∑_i∑_ij⟨⟩ł()̊ϵł()̊
                
                Raman Research Institute, Bangalore 560080, India
                
                
                Centre for Condensed Matter Theory, Department of Physics,
                 Indian Institute of Science, Bangalore 560012, India
                
                madan@ncbs.res.in
                on lien from : Raman Research Institute, C.V. Raman Avenue, Bangalore 560080, India
                
                
                Simons Centre for the Study of Living Machines, National Centre for Biological Sciences (TIFR), Bellary Road, Bangalore 560065, India
                
                Recent experimental studies, both  in vivo and  in vitro, have revealed that membrane components that bind to the cortical actomyosin meshwork are driven by active fluctuations, whereas 
                membrane components that do not bind to cortical actin are not.
                 Here 
                we study the statistics of density fluctuations and dynamics of particles advected in an active  quasi-two dimensional medium 
                comprising self-propelled filaments with no net orientational order, using a combination of agent-based Brownian dynamics simulations and analytical calculations. The particles interact with each other and with the
                self-propelled active filaments via steric interactions. 
                We find that the particles show a tendency to cluster and their density fluctuations reflect their binding to and driving by the active filaments.
                The late-time dynamics of tagged particles 
                is diffusive, with an  active diffusion coefficient that is independent of (or at most weakly-dependent on) temperature 
                at low temperatures. Our results are in qualitative agreement with 
                the experiments mentioned above. In addition, we make predictions that can be tested in future experiments.
                47.63.mh, 87.17.-d, 05.60.cdInteracting passive advective scalars in an active medium
                
                    Madan Rao
                    
                ==========================================================
                
                
                
                § INTRODUCTION AND MOTIVATION
                
                
                
                The spatial organization, clustering and dynamics of many cell surface molecules is influenced by interaction with 
                the actomyosin cortex <cit.>, a thin layer of actin cytoskeleton and myosin motors, measured to be around 250 nm in HeLa cells <cit.>, that is 
                juxtaposed with the membrane bilayer. Although the ultrastructure of the cortical actin cytoskeleton is as yet poorly defined, there is growing
                evidence that it is composed simultaneously of dynamic filaments <cit.> and an extensively branched static  meshwork <cit.>. 
                The coupling of the membrane to these two types of actin configurations is expected to affect the dynamics and
                organization of membrane components, presumably in different ways.
                
                A common feature of these cell surface proteins is that they can bind, directly or indirectly, 
                to cortical actin. As a consequence, the action of myosin motors on actin at the cortex help drive the local clustering and dynamics of these cell surface proteins.
                Mutations of these proteins that abrogate this actin binding capacity, leave them unaffected 
                by the dynamics of the actomyosin cortex <cit.>. 
                Similarly, silencing myosin motor activity, renders the dynamics of these cell surface proteins normal and akin to their mutated counterpart <cit.>.
                 A description of the cell surface as an  Active Composite of
                a multicomponent, asymmetric bilayer juxtaposed with a thin cortical actomyosin layer (Fig. <ref>), appears to
                consistently explain the anomalous dynamical features of these proteins <cit.>. In these earlier studies 
                 <cit.>, we had used a coarse grained description, based on active hydrodynamics <cit.>.
                Recently, it has been shown that much of this behaviour is recapitulated in a minimal 
                 in vitro reconstitution of a supported bilayer in contact with a thin layer of short actin filaments 
                and Myosin-II minifilaments, driven by the hydrolysis of ATP <cit.>. The success of this minimal setup motivates us 
                to revisit our continuum hydrodynamic description from a more microscopic standpoint.
                Our present study is an agent-based Brownian dynamics simulation of a mixture of
                polar active filaments and passive particles which interact with each other. Since the parameter space for exploration is large, 
                we restrict our study here to the case where the density of filaments and particles is low (negligible filament overlap), 
                and where there is no net orientational order.
                
                
                For the present purposes, we schematically represent the Active Composite Cell Surface, as in Fig. <ref>.
                To enable a systematic study, we conceptually separate out the different architectures of actin at the cortex -
                (i) where the cortical layer consists of only short dynamic filaments described as an active fluid (the subject of the 
                present paper) and (ii) where the cortical  layer consists of only long filaments forming a static mesh of 
                characteristic mesh size ξ≈ 52 nm in FRSK cells <cit.> (which we take up in a later publication). 
                In future, we will combine these configurations into a single model of the cortical actomyosin.
                
                
                
                
                 Our choice of Brownian dynamics simulations
                  is motivated by experiments on tagged particle diffusion both on the cell surface and in the  in vitro reconstitution.
                 Molecules that bind to dynamic actin (passive molecules) are affected by the 
                active fluctuations of actomyosin - their diffusion shows anomalous behaviour strongly indicative of active driving. On the other hand, 
                molecules that do not interact with actin (inert molecules), such as short chain lipids and proteins whose actin-binding domain has been 
                mutated so as to abrogate their interaction with actin, do not show any influence  of active fluctuations <cit.>. There appears to be no sign that the transport of these
                inert molecules is affected by potential 
                hydrodynamic flows induced by active stresses coming from actomyosin <cit.>. 
                
                
                While our primary motivation are the experimental studies of the tagged particle dynamics on the cell surface, our work is also relevant to transport in other
                living and nonliving systems, as long as the effects of hydrodynamics are negligible, for instance, to the movement of multiple motor-driven cargo vesicles or synthetic beads on the cytoskeletal 
                network <cit.>.
                
                
                
                
                
                § BROWNIAN DYNAMICS AND CHARACTERIZATION
                
                
                 §.§ Simulation details
                
                
                We study the dynamics of a mixture of polar active filaments and passive particles using 
                an agent-based Brownian dynamics simulation. 
                The passive particles are modelled as mono-disperse soft spheres of diameter σ. A pair of passive particles separated by a distance r interact via a truncated Lennard Jones (LJ) pair potential
                of the form,
                
                V_pp(r)     =     4 ϵ[(σ/r)^12-(σ/r)^6] + V_0 +V_2 r^2          r≤ r_c  
                    =        0            r > r_c 
                        
                where r_c=2.5 σ and values of V_0 and V_2 are chosen so that the potential and force are continuous at the truncation point. 
                We set σ=1 and ϵ=1 to be the units of length and energy, respectively.
                
                The polar filaments are modelled as semi-flexible bead-spring polymers, with 
                both stretch and bend distortions.
                We implement excluded volume interaction between the beads of same filament,
                as well as between two different filaments through a truncated Lennard Jones pair potential of the form,
                
                V_bb(r)     =     4 ϵ^'[(σ^'/r)^12-(σ^'/r)^6] +V_0^'+ V_2^' r^2        r≤ r^'_c  
                    =     0           r > r^'_c
                           
                where 
                r is the distance between the centres of the corresponding beads, 
                r^'_c=2^1/6σ^' and V_0^', V_2^' are constants, 
                chosen so that the potential and force are continuous at r^'_c.
                We take σ^'=2 and ϵ^'=1. Each filament is composed of
                10 beads and therefore has an equilibrium length l_0=20, in the units of σ.
                
                Note that with our choice of cutoffs, the particle-particle interaction V_pp has both attractive and repulsive parts, whereas the 
                bead-bead interaction V_bb is strictly repulsive.
                
                To make the filament semi-flexible, we impose additional spring forces on the beads. 
                A harmonic stretching potential with 
                spring constant K_c=400, in units of ϵ/σ^2, ensures that the length of the filament does not deviate significantly from its
                equilibrium value, l_0=20. The bending energy of a triplet
                of connected beads is also harmonic in the angle, with a bending stiffness 
                K_a = 600, in units of ϵ. This high K_a makes the filaments very stiff, with 
                a typical persistence length much larger than l_0.
                
                A propulsion force f=f_0 n̂ is imposed 
                on each of the beads, along the average direction (n̂) of all the bonds present in a filament. 
                Note that we do not impose any filament alignment rule nor do we prescribe any activity decorrelation time. 
                Instead, these originate from thermal fluctuations on the constituent monomers comprising each filament and 
                collisions driven by thermal and active forces, an emergent many-particle feature. As a consequence, 
                both the local alignment and orientational de-correlation time are functions of 
                temperature, density  and activity.
                
                The interactions between the beads of the filament and the passive particles are modelled by a 
                harmonic potential of spring constant K_s=50 in units of ϵ/σ^2. The harmonic potential is
                truncated at a cutoff distance r_0= 1 and set to zero beyond it. When a passive
                particle comes within a distance r_0 from the centre of a filament bead, it binds 
                to the corresponding bead and gets advected along with the filament, under the application of propulsion force f_0 n̂.
                
                
                Unless mentioned otherwise, all results presented here are for N_p=800 passive particles 
                and N_r=50 self-propelled filaments in a two dimensional area of linear dimension L=396.4. For most of the study, we
                take the area fractions of the rods (c) and particles (ρ) to be c = 0.01 and ρ=0.004, respectively.
                
                
                
                The Brownian dynamics equations involve an update of both the passive particle and the filament bead coordinates, for which we have used a simple Euler integration scheme with
                integration time step Δ t ∼ 10^-4.
                The dynamics of the particle coordinates is given by
                
                
                ṙ_i^p=
                    
                      -γ^-1_p ∇_i V_p   +  √(2 k_BT/γ_p)ξ_i,    (unbound) 
                
                      -γ^-1_p ∇_i V_p   +  √(2 k_BT/γ_p)ξ_i+  f_i/γ_p,    (bound)
                where
                γ_p is the friction coefficient of the passive particle, V_p is the net potential felt by the i-th passive particle and includes contributions 
                from Eq. <ref> and the bead-particle spring interactions.
                The diffusion of the unbound particle is driven by a thermal noise
                ξ_i  with zero mean and unit variance acting on i-th particle.
                On the other hand, the bound particles are subject to both the thermal noise and active driving.
                 
                The dynamics of the filament-bead coordinates is given by
                ṙ_j^b = -γ^-1_b ∇_j V_b   + √(2 k_BT/γ_b) ξ_j +   f_j/γ_b                                         
                where
                γ_b is the friction coefficient of the bead, V_b is the net potential felt by the j-th bead and includes contributions from Eq. <ref>, harmonic stretching and bending interactions.
                
                
                 We take γ_p=1, which together with σ=1 and ϵ=1, sets the units of space, time and energy. 
                 All other quantities can be written in terms of these units, so as to make Eqs. <ref>, <ref> dimensionless.
                In all that follows below,  except in Sec. V,  we have taken γ_b=γ_p.
                 To be able to make contact with experiments, we have to translate our simulation units (S.U.) to
                real units (R.U.). Setting σ = 10 nm,  γ=0.123 pN μ m^-1 s <cit.> and ϵ = 4.14 × 10^-21 J, 
                we can convert our simulation units to real units, as displayed in Table 1.
                 
                We have typically run the Brownian dynamics simulation for a total time t ∼ 10^4, ensuring that 
                the system has reached steady state.
                We have varied the temperature over T = 0.25, 0.5, 0.75, 1, 2, 4, 6, 8, 10, 15, 25 and the
                propulsion force over f_0 = 0.0, 0.5, 1.0, 2.0, 3.0, 4.0. Our initial conditions are chosen
                from a thermal distribution at temperature T and all results presented here are averaged over 16 such independent
                initial realisations.
                
                
                
                
                 §.§ Statistics of filament orientation
                
                We characterise the i-th filament by its  centre of mass position r_i and a unit vector n_i =(cosθ_i , sinθ_i ) 
                along its long axis to describe its polar orientation (recall that the filaments are very stiff). 
                We first ensure that the configuration of filaments is in the spatially homogeneous, orientationally isotropic state - this is demonstrated in 
                the plots of the probability distribution of the polar P(θ) and nematic orientations P(θ̃) (Fig. <ref>).
                
                
                
                We then calculate the orientational correlation lengths, so as to ensure that this is much smaller than our system size and comparable to the size of the filaments. To do this, we
                 calculate the spatial correlations  of both the polar and nematic orientation,
                
                C_P(r)     =    ⟨1/N^2∑_i=1^N ∑_j=1^N  cos(θ_i -θ_j) ⟩
                
                C_N(r)     =    ⟨1/N^2∑_i=1^N ∑_j=1^N ( 2  cos^2(θ_i -θ_j) -1 ) ⟩
                where r= | r_i -  r_j | is the distance between the centre-of-mass of the i-th and j-th filaments.
                By fitting this to an exponential (Figs. <ref> and <ref>), we extract the polar and nematic orientation correlation lengths, ζ_P and ζ_N,
                whose dependence on the area fraction of  filaments  c (we will henceforth refer to this as filament density) is shown in the inset.
                
                
                
                
                Throughout the paper (unless mentioned otherwise) we work at a filament density of c=0.01, and temperatures T≥0.5; in this regime, the
                orientational correlation lengths are of the order of the filament length, l_0, and hence comfortably within the isotropic phase.
                
                
                
                 §.§ Statistics of filaments persistence
                
                A self-propelled filament moves persistently in a direction set by V_a=f_0 n̂/γ_b, taken to be along its polar orientation, until its orientation gets decorrelated, either due to thermal noise on the individual monomers constituting the
                filaments or due to inter-filament 
                collisions. What is the spatial correlation of these directions of persistent motion? Here, we calculate this velocity orientation correlation function C_V(r) at different T (Fig. <ref>) using the formula
                
                C_V(r)     =    ⟨1/N^2∑_i=1^N ∑_j=1^N ( 2  cos^2(Θ_i -Θ_j) -1 ) ⟩
                where Θ_i represents the direction of the velocity vector of the i-th filament.
                From this spatial correlation function C_V(r) we extract a 
                correlation length (ζ_V) which measures the spatial extent of the persistent dynamics. We find again that at c=0.01 over the temperature range T ≥ 0.5,
                 the correlation length ζ_V is comparable with the filament length l_0.
                 
                
                
                 §.§ Statistics of binding-unbinding of passive particles onto filaments
                
                
                The dynamical equations (<ref>) and (<ref>) are written entirely in terms forces, either active or derived from a potential, and thermal noise. The particles experience a 
                binding and unbinding onto the filaments which depend on this interplay between thermal noise and the attractive potentials. Thus for instance, 
                the unbound passive particles diffuse in the two dimensional medium and ever so often come within the vicinity (r ≤ r_0 = 1) of a moving filament-bead, whereupon they bind to the filament-bead.
                In the low density limit, we expect the binding rate k_b to be diffusion limited and so k_b ∝ T and independent of K_s where K_s is the strength of trapping harmonic interaction. 
                
                
                To study the unbinding of a particle bound from a filament-bead, we compute the rate of escape of a particle trapped in a truncated attractive harmonic potential <cit.>,
                parametrised by K_s and r_0. This is given by
                
                k_u = K_s^2 r_0^2/γ_p k_B T exp(-  K_s r_0^2/2 k_B T)
                
                and should be a good description of the dynamics of unbinding of the passive particles in the limit of low particle density.
                
                
                We compare these theoretical estimates with the results of simulations
                on a mixture of particles and filaments at equilibrium (no active propulsion), from which
                we extract the values of k_u and k_b in two different ways. In the first method,
                we represent the stochastic binding and unbinding by a 
                telegraphic process <cit.>, characterised by a mean duty ratio,
                ϕ = k_b /k_b + k_u  ,
                
                the fraction of
                time spent by the tagged particle in the  bound state over the observation time, and a 
                 two-point correlator,
                ϕ(t) ϕ(t^')  =  ϕ^2 +  ϕł 1- ϕe̊^-2 | t-t^'|/t_sw
                where 
                 t_sw = 2/k_b+k_u  ,
                
                is called the mean switching time and describes the mean time taken to switch from a bound to an unbound state.
                We  calculate k_b and k_u, by fitting our simulation results to  ϕ  and ϕ(t) ϕ(t^'). In the second method,
                we calculate k_u (k_b) directly, from the inverse mean time that the particle stays bound (unbound) on the filament.
                
                The two numerical methods give identical results (after scaling by a constant factor) as a function of the particle-filament  binding potential K_s and T, and that these
                agree with the analytical estimates, Fig. <ref>, with no fit parameter.
                
                It is important that we do not  prescribe the binding-unbinding rates, rather 
                we derive 
                them from the assigned potentials. The binding-unbinding rates thus depend nontrivially on temperature;  they would also depend on the
                density of passive particles and filaments in the high density limit.
                This will  be crucial to our estimation of the tagged particle diffusion coefficient and its comparison with experimental data (Sect. III D).
                
                
                
                § DENSITY FLUCTUATIONS OF PASSIVE ADVECTIVE SCALARS IN AN ACTIVE MEDIUM
                
                
                We now study the statistics of density fluctuations and dynamics
                of the actively driven passive particles. We find that the active driving tends to cluster the passive particles; this shows up in
                 the two point spatial density correlation function and the statistics of the density fluctuations. 
                
                
                
                 §.§ Radial distribution function
                 
                
                
                We study the behaviour of the radial distribution function of the passive particles g(r),
                
                    g(r)= 1/N_pρ⟨∑^_i∑^_i≠ jδ(r-| r_i - r_j |) ⟩,
                
                where N_p is the total number of passive particles and ρ  the passive particle density. 
                When k_b=0 and f_0=0, i.e., when the particles do not bind to the filament (inert particles) and there is no propulsion force, g(r) has the form of a dilute fluid (Fig. <ref>).
                When we allow for particle binding, but in the absence of propulsion force, the g(r) displays 
                oscillations, which arise from particles binding to periodic locations on the filaments (Fig. <ref>) - note r=20 coincides with the filament  length. In this equilibrium situation, the 
                 particles not bound to the filaments do not show any clustering.
                
                We now consider the case when the filaments are driven by a propulsion force f_0. We find that in addition to the periodically spaced bound particles, there is a significant fraction of unbound particles that are clustered
                 (Fig. <ref>).
                This clustering is a consequence of the shepherding of  passive particles by the active filaments - this has also been reported in <cit.>, with the difference that in that case the dynamics of
                the particles affects the active filaments.
                
                
                
                 §.§ Probability distribution of local number density
                 
                This activity induced clustering of the passive particles should be reflected in the probability distribution of the excess number density.
                 To  compute this we
                divide the  system into blocks of size Ω = 39.64 and
                count the number of passive particles n in each block, to obtain the 
                 steady state distribution P(n).
                
                
                The inert particles, in the dilute limit (Fig. <ref>), has a probability distribution that resembles a
                gas at temperature T and the average number of particle in the blocks n̅, namely, 
                
                    P(n)=1/√(2 πσ^2)exp[ -(n- n̅)^2/2 σ^2]
                 
                where n̅=Ωρ. We have taken the variance of the above distribution to have a virial form, σ^2 = n̅ - 2 B_2/Ωn̅^2 + … - a fit to the numerical data
                gives B_2=-0.99 (dark line in Fig. <ref>).
                
                
                
                
                On the other hand, the probability distribution for passive particles picks up an exponential tail arising from the binding-unbinding statistics of the particles.
                This exponential tail gets more pronounced when the filaments are made active, and 
                which moves towards the typical value as the active propulsion force gets larger (Fig. <ref>). 
                This reflects the fact that for high driving, the typical particle is clustered. Both these results are entirely consistent with 
                the experimental observations reported in Fig. 4 D of <cit.>.
                
                
                
                
                 §.§ Number fluctuations : crossover from anomalous to Brownian
                 
                Note that the active system of filaments is in the isotropic phase and we should not expect to see giant number fluctuations normally associated with
                active systems with global orientational order <cit.>. However when we compute the root mean square fluctuations △ n and mean 
                n̅ of the number of passive particles over regions of ever increasing area, and plot them with respect to each other, we find that initially △ n ∝n̅^α
                with α = 0.784.
                 Subsequently, as n̅ increases, the variance scaling shows a cross over to α=0.5.
                This crossover occurs over a scale corresponding to the orientational correlation length, which can in principle be large, especially close to the
                isotropic-nematic transition or high f_0. 
                This is especially apparent in the high particle density regime, see Fig. <ref> for particle density ρ=0.05 and filament density c=0.02.
                This slow crossover is
                 consistent with the experiments on the  in vitro reconstitutions of actomyosin on a supported bilayer <cit.>.
                
                
                
                § TRANSPORT OF PASSIVE ADVECTIVE SCALARS IN AN ACTIVE MEDIUM
                
                
                We now study the transport of passive particles moving in the active medium. Because the filaments are orientationally disordered, the long time dynamics of the particles
                is always diffusive. However the diffusion characteristics can change depending on the statistics of (un)binding to the active filaments.
                
                
                
                 §.§ Typical trajectories
                
                The space-time trajectories of the passive particles show
                three qualitatively different behaviours. At very low temperatures compared to U = K_sr_0^2/4, 
                a passive particle once bound
                to a filament, rarely unbinds, and hence gets advected with the self-propelled filament (Fig. <ref>(a)).
                
                The direction of 
                advection changes because of thermal fluctuations and collisions between filaments. 
                
                Increasing the temperature increases the probability of unbinding from the filament, whereupon the particle undergoes unrestricted thermal diffusion before 
                binding again (Fig. <ref>(b)). At even higher temperatures, k_BT/U≈  1  the particles do not bind to the filaments and the motion is simple thermal diffusion
                (Fig. <ref>(c)).
                
                
                
                 §.§ Statistics of displacements and diffusion coefficient
                 Propensity distribution.
                The distribution of displacements  Δ x (along the  direction) evaluated over a time window t_w is called the
                 propensity distribution. This will depend on the statistics of binding/unbinding, which
                in turn depends on the temperature T and f_0, densities of filaments and particles, and of course on the time window t_w, which we fix at t_w=10. This can be obtained both from our Brownian dynamics 
                simulation and, in the dilute limit, analytically.
                
                In the dilute limit, one can obtain the form of this probability distribution from the stationary process describing  the particle vector-displacements in a small time interval t_w,
                	 r(t_w)  =   r(0) + ∫_0^t_w(t')  dt' 
                	
                	where   is the velocity of the tagged particle at time t, given by,
                		(t)  =  ϕ(t)  f_0/γ n̂(t) + (1-ϕ(t))  ξ(t),
                
                where n̂ is the polar vector representing the orientation of the filament to which the particle is bound at time
                t, ξ is the thermal noise, and 
                ϕ(t) is the  telegraphic noise whose statistics is described in Sect. II D. The distribution of the particle displacements Δ x_t_w≡ ( r(t_w) -   r(0) )·,
                can be obtained by evaluating,
                
                P(Δ x) = ⟨δ(Δ x-Δ x_t_w) ⟩
                where Δ x_t_w is obtained from Eq. <ref> and the angular bracket denotes an average over the joint distribution of ξ and ϕ. This can be evaluated by standard
                techniques of Fourier transformation and cumulant expansion <cit.>,
                lnP̃(k) = ∑_m=1^∞(i k)^m/m !⟨(Δ x_t_w)^m  ⟩_c
                
                where  P̃(k) is the Fourier transform of P(Δ x). The m-th cumulants ⟨(Δ x_t_w)^m  ⟩_c can be evaluated from Eqs. <ref>, <ref>, knowing 
                 that ϕ and ξ are independent stochastic processes. The distribution P(Δ x) is then obtained by taking the inverse Fourier transform of  P̃(k).
                
                However, in practice, the inverse Fourier transform of P̃(k), for the stationary process Eq. <ref>, has to be evaluated numerically.
                Rather than do this, we 
                 provide an alternate argument which gives more insight.
                
                At low enough T, the 
                passive particles are completely bound to the self-propelled filaments, and so as long as t_w < τ, the 
                orientational
                correlation time of the filaments,
                the particles get displaced by Δr=f_0 t_w/γ_pn̂, where 
                γ_p is the friction coefficient and n̂ is the unit vector representing the average orientation of a filament during time interval
                t_w. 
                Since the filament orientation is uniformly distributed, the contribution to the step-size distribution from this process is 
                P(Δ x)=γ_p/π f_0 t_w1/√(1-(γ_p/f_0Δ x/t_w)^2). On the other hand, at very high T, the particles are completely unbound and undergo thermal diffusion,
                for which the step-size distribution is P(Δ x)=1/√(2 πσ^2)exp[-
                (Δ x)^2/2 σ^2], where σ^2=2  k_B T/γ_p t_w. 
                We propose that at an intermediate T, the propensity distribution can be written as a  linear combination, weighted  by 
                 the duty ratio ⟨ϕ⟩, i.e., 
                
                P(Δ x)= γ_p/π f_0 t_w⟨ϕ⟩/√(1-(γ_p/f_0Δ x/t_w)^2) + 
                1-⟨ϕ⟩/√(2 πσ^2)exp[-(Δ x)^2/2 σ^2].
                
                The agreement of this approximate analytical form, with no undetermined parameters, with the results of the Brownian simulation is quite reasonable, see Fig. <ref>.
                
                 Mean square displacement. From the statistics of the displacement we can compute the 
                 mean square displacement (MSD) as ⟨Δ r^2(t) ⟩=⟨1/N_p∑_i | r_i(t_0+t) -  r_i(t_0)|^2  ⟩, where r_i is the position of i-th particle. This shows a change from a short time diffusive regime crossing over to a long time diffusive regime via an
                 intermediate super-diffusive regime (Fig. <ref>(a)). We estimate the second crossover time t_c(T,f_0) from the super-diffusive to late time diffusion D, by fitting the simulation data to ⟨Δ r^2(t) ⟩= 4Dt [1-exp(-t/t_c)] <cit.>,
                using which
                we can collapse the MSD data for different values of active propulsion f_0 (Fig. <ref>(b)). From this we see that t_c decreases with
                f_0 (Fig. <ref>(b) inset). This is because the filament orientation decorrelates on account of collisions, whose frequency increases with f_0.  In experimental systems  where  hydrodynamics plays a crucial 
                role <cit.>, this dependence of t_c on
                f_0 may be different.
                
                In situations where the crossover t_c is large, the apparent super-diffusion behaviour would last for many decades in time.
                We can then fit the MSD to  ⟨Δ r^2(t) ⟩∼ t^β to obtain an super-diffusion exponent β > 1 - 
                we find that β=1.95 at T=0.5 and β=1.77 at T=2.0, for a propulsion force f_0=4.0.
                 Temperature and activity dependence of MSD.
                We now compute the late time diffusion coefficient of the tagged particles, D = lim_t →∞⟨Δ r^2(t) ⟩/4t, for different T and f_0. 
                For a fixed f_0, one might expect that at low temperatures D is weakly dependent on (or even independent of)
                temperature because a particle once bound to the filament remains so and undergoes  active diffusion as it is transported by the filament
                (Fig. <ref>). 
                As we increase the temperature, D decreases, since a particle spends less time, on an average, bound to the filament (recall we have set γ_b=γ_p).
                At high temperatures, the particles are predominantly unbound, and hence  D resembles that of an inert particle, which increases linearly with temperature.
                This is indeed what we see from a direct numerical simulation of the Brownian dynamics trajectories of a tagged particle (Fig. <ref>).
                
                
                
                
                From the stationary process, Eq. (<ref>), the MSD of the tagged passive particle,
                δ r^2(t)   =   ∫_0^t ∫_0^t  (t') ·(t”) dt' dt”  ,
                
                 immediately gives the diffusion coefficient, 
                
                D = 1/2∫^∞_0(t) ·(0) dt.
                
                
                Using Eq. (<ref>), we see that the diffusion coefficient of the bound particle is given 
                by the correlations of _a, which  is given by (Fig. <ref>(a)),
                 _a(t) ·_a(0)     =     f_0^2/γ_pcos(θ(t)-θ(0)) 
                    =     f_0^2/γ_p e^- | t|/τ
                  The diffusion coefficient can now be simply evaluated,
                  
                  D=f_0^2 τ k_b(τ k_b +1)/2 γ_p^2(k_u+k_b)[(k_u+k_b)τ+1]+k_B T k_u/γ_p(k_u+k_b)
                
                  
                  To plot D versus T and f_0, we need to know the values of k_b, k_u (equivalently ⟨ϕ⟩, t_sw) and τ, which depend on the temperature and density, and which
                  we obtain from our simulations.
                  We then compare this semi-analytical form to the direct numerical computation of the diffusion coefficient from the Brownian dynamics trajectories (Fig. <ref>).
                  The agreement between the two is excellent. These observations are entirely consistent with the results of <cit.>.
                  
                  
                
                
                  
                  It might be objected that in our analysis we have treated the active propulsion as an independently tunable parameter, thus precluding the possibility that the 
                  activity itself may be temperature dependent. However, as we saw in <cit.>, and as noted elsewhere <cit.>, the actomyosin contractile processes taken as a whole, 
                  appear to be independent of temperature in the physiological range,
                24^∘-37^∘C.
                
                Figure <ref> shows the dependence of D on filament concentration c, both from direct simulations and from the analytical form using the values of 
                k_b, k_u (equivalently ⟨ϕ⟩, t_sw) and τ,
                from simulations. This shows optimal transport at a specific filament concentration; the orientational decorrelation time is smaller at higher filament concentration, due to higher collision
                frequency.
                
                
                
                § VISCOSITY STRATIFICATION AND ITS EFFECT ON MEMBRANE DIFFUSION
                
                So far, our study of transport of passive molecules in an active medium has been restricted to two dimensions. However as we discussed in Sect. I, the cell surface is a composite of a bilayer membrane and a 
                thin actomyosin cortex. Thus while the proteins move on the cell membrane, the actively driven actin moves in the actomyosin cortex. The viscosities of these two layers are significantly different, with the
                bilayer membrane having a viscosity which is an order of magnitude larger than the cortex (≈0.86 Pa s <cit.>).
                Indeed the local viscosity of a multicomponent membrane can be quite heterogeneous - for instance the particle mobility within the so-called “membrane rafts” or liquid-ordered regions on the cell membrane
                can be very different from those within liquid-disordered regions. Moreover the local cortical viscosity depends on local actin, myosin and cross-linker concentrations.
                 How does this viscosity mismatch affect the actively driven transport of passive molecules?
                
                To address this important issue within our simulation, we vary the ratio of the friction coefficients Γ = γ_p/γ_b in Eqs. <ref>, <ref>. 
                We find that the mean fraction of passive particles bound to filaments ⟨ n_b ⟩
                 decreases with increasing Γ over a range of T and f_0=4.0 (Fig. <ref>(a)).
                This is an interesting observation, since one might have naively thought that ⟨ n_b ⟩ is solely governed by binding-unbinding, a purely equilibrium process and hence 
                independent of relative viscosities. However, we see that the drag induced by the imposed viscosity stratification (a nonequilibrium feature), can “peel-off” particles from the filaments.
                It is not clear to us why we see a shoulder at intermediate values of Γ for low enough temperatures (Fig. <ref>(a)).
                
                
                
                This is reflected in changes that we observe in the measured diffusion coefficient D, as it decreases with increasing Γ at different T (Fig. <ref>(b)).
                As can be seen, the active-diffusion regime at low temperatures becomes significantly more temperature dependent as the viscosity mismatch Γ increases.
                
                The results of this section are not purely academic, on the contrary taken together they pose an interesting possibility that by tuning local viscosity mismatch, 
                for instance by locally recruiting the so-called “membrane rafts” or liquid-ordered regions on the cell membrane or by  locally regulating the concentrations of actin, myosin or cross-linkers,
                the living cell surface could control  the clustering and transport of
                specific membrane proteins.
                
                
                
                § DISCUSSION
                
                We had earlier shown that a coarse grained active hydrodynamics description of the active composite cell surface, successfully explains the statistics of clustering 
                of membrane proteins capable of binding to the cortical actomyosin in living cells <cit.>. Such a description make predictions regarding the statistics of density fluctuations and
                transport of such actin-binding membrane proteins, which were verified in experiments <cit.>. Following this we were able to recapitulate much of this behaviour in a minimal
                 in vitro system comprising  a thin layer of short actin filaments and Myosin-II minifilaments on a supported bilayer <cit.>. The success of this approach has motivated us to do a
                agent-based Brownian dynamics simulation using these minimal ingredients - that of a collection of passive molecules which bind/unbind to actin filaments and move in this active medium in two dimensions.
                
                
                The  results obtained  here, based on simulations and analytical calculations,
                 are in qualitative agreement with 
                the experiments both  in vivo and  in vitro. For instance, the exponential tails appearing in the probability distribution of the number (Fig. <ref>) and the 
                scaling of the variance of the number (Fig. <ref>) is precisely the behaviour seen in our earlier  in vitro experiments.
                In addition, we show how activity induced clustering of passive particles (Fig. <ref>) arises naturally from such a minimal description.
                
                We have also studied transport of passive particles moving in this active medium, and find that there is a crossover from an intermediate time super diffusive to late time
                diffusive behaviour as a consequence of active driving (Fig. <ref>(a)). The transport behaviour shows a striking  dependence on temperature and active forcing - at low temperatures
                the diffusion coefficient is insensitive to temperature, and crosses over to a linear temperature dependence at higher temperatures, in qualitative agreement with 
                experiments <cit.>.
                
                
                Finally, recognising that the viscosity of the cortical layer is different from that of the membrane, we show that a friction coefficient mismatch has a strong effect on the mean
                number of bound particles and the diffusion coefficient. This is a consequence of the drag induced by the imposed viscosity stratification, which results in a 
                 “peeling-off” of the particles from the filaments. This opens up the possibility of  local tuning of viscosity mismatch, for instance by locally recruiting the so-called “membrane rafts” or
                  liquid-ordered regions on the cell membrane or by locally regulating the concentrations of actin, myosin or cross-linkers. This could result in yet another mechanism by the cell surface might locally
                 control the clustering and transport of specific membrane proteins. We hope that some of these predictions can be tested in future experiments.
                
                
                
                
                § ACKNOWLEDGEMENT
                
                We thank K. Husain, R. Morris and D. Banerjee for critical inputs. 
                MR thanks S. Mayor and D. Koster for long years of collaboration.
                SKRH and RM thank the Simons Centre at NCBS for computational facilities and hospitality.
                RM acknowledges financial support from CSIR, India.
                
                
                raomayor  M. Rao and S. Mayor, Curr. Op. Cell Biol.  29, 126 (2014).
                  
                  Goswami D. Goswami et al., Cell  135, 18085 (2008).
                  
                  kripa K. Gowrishankar, S. Ghosh, S. Saha, S. Mayor and M. Rao, Cell  149, 1353 (2012).
                  
                  saha S. Saha, I. Lee, A. Polley, J. T. Groves, M. Rao and S. Mayor, Mol. Biol. Cell  26, 4033 (2015).   
                  
                  Hancock S.J. Plowman, C. Muncke, R.G. Parton and J.F. Hancock, Proc. Natl. Acad. Sci. USA  102, 15500 (2005).
                  
                  Lingwood D. Lingwood and K. Simons, Science  327, 46 (2010).
                  
                  Paluch P. Chugh, A.G. Clark, M. B. Smith, D. A. D. Cassani, K. Dierkes, A. Ragab, P. P. Roux, G. Charras, G. Salbreux and E. K. Paluch, Nat. Cell Biol.  19, 697 (2017).
                                    
                  Morone Morone et al., J. Cell Biol.,  174, 851 (2006).
                  
                  kripasoftK. Gowrishankar and M. Rao, Soft Matter  12, 2040-2046 (2016).
                
                  rmpM. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys.  85, 1143 (2013).
                  
                   dariusD.V. Koster, K. Husain, E. Iljazi, A. Bhat, P. Bieling, R.D. Mullins, M. Rao and S. Mayor, Proc. Natl. Acad. Sci. USA  113, E1645 (2016).
                
                  fujiwara T. Fujiwara, K. Ritchie, H. Murakoshi, K. Jacobson and A. Kusumi, J. Cell Biol.  157, 904 (2002).
                  
                    yhatY. Hatwalne, S. Ramaswamy, M. Rao, and R. A. Simha, Phys. Rev. Lett.  92, 118101 (2004).
                  
                  abasuA. Basu, J. F. Joanny, F. Julicher, and J. Prost, New J. Phys. 14, 115001 (2012).
                  
                  GrannickK. Chen, B. Wang, and S. Granick, Nat. Mater.  14, 589 (2015).
                  
                    kamm T. C. Bidone, W. Jung, D. Maruri, C. Borau, R. D. Kamm and T. Kim, PLoS Comput Biol  13, 1005277 (2017).
                  
                
                  
                  grebenkovD. Grebenkov, J. Phys. A 48, 013001 (2015).
                  
                  kampenN.G. Van Kampen in  Stochastic Processes in Physics and Chemistry, 3rd Edition, North-Holland, Amsterdam (2007).
                  
                  nitin N. Kumar, H. Soni, S. Ramaswamy and A. K. Sood, Nat. Comm.  5, 4688 (2014).
                  
                  toner-tu J. Toner and Y. Tu, Phys. Rev. Lett.  75, 4326 (1995).
                  
                  vijay V. Narayan, S. Ramaswamy and N. Menon, Science  317, 105 (2007).
                    
                  simhaR. A. Simha, J. Toner and S. Ramaswamy, EPL  62, 2 (2003).
                    
                
                  
                
                  
                  libchaber X. Wu and A. Libchaber, Phys. Rev. Lett.  84, 3017 (2000).
                
                
                  
                  spudich M. P. Sheetz, R. Chasan and J. A. Spudich, J. Cell Biol. 99, 1867-1871 (1984).
                  
                  levy H.M. Levy, N. Sharon, and D. E. Koshland Jr.,  Proc. Natl. Acad. Sci. USA  45, 785-791 (1959).
