Entropic Effects of Interacting Particles Diffusing on Spherical Surfaces

We present a molecular dynamics and theoretical study on the diffusion of interacting particles embedded on the surface of a sphere. By proposing five different interaction potentials among particles, we perform molecular dynamics simulations and calculate the mean square displacement (MSD) of tracer particles under a crowded regime of high surface density. Results for all the potentials show four different behaviors passing from ballistic and transitory at very short times, to sub-diffusive and saturation behaviors at intermediary and long times. Making use of irreversible thermodynamics theory, we also model the last two stages showing that the crowding induces a sub-diffusion process similar to that caused by particles trapped in cages, and that the saturation of the MSD is due to the existence of an entropic potential that limits the number of accessible states to the particles. By discussing the convenience of projecting the motions of the particles over a plane of observation, consistent with experimental capabilities, we compare the predictions of our theoretical model with the simulations showing that these stages are remarkably well described in qualitative and quantitative terms.


INTRODUCTION
In several physical and biological systems the mass transport phenomena is carried out in surfaces with non vanishing curvatures. The diffusion of bio-molecules and other particles on the surface of liposomes, drops and other curved entities is of high relevance due to its potential applications in biomedicine and technology [1][2][3][4]. Characteristic examples of this phenomena can be found in the diffusion of proteins on curved membranes or in the surface diffusion of molecules and chemical solvents in catalytic surfaces [5][6][7][8].
In these curved domains, many geometric aspects relating the motion of the particles on surfaces have been studied in recent years [9][10][11][12][13][14]. These studies are relevant in the understanding of applied problems such as nucleation, spinodal decomposition, adsorption, and phase transitions, since the dynamics of the particles is modified due to the curvature of the surface on which the molecules are embedded [15][16][17][18].
In particular, the most relevant quantity to measure in these particle motions is the mean square displacement (MSD). However, most of the previously refereed studies provide geometric information of the MSD and the effective diffusion coefficient only for the case where the interaction of the particles occurs with the surface and not among particles themselves. In this simplification, many expressions can be found in the literature for approximating the short and long time behavior of the MSD, see Ref. [9] and references therein. However, for the case of interacting particles, there are only few numerical works studying the influence of interaction and confinement in the MSD [19][20][21].
In this work we will consider how the surface diffusion displacements are influenced by the interaction of the particles and the projection of their displacement in a plane of observation. This is because for most experimental cases of interest, the geodesic trajectories of the particles cannot be followed at all times, and only the projection in a plane of observation can be measured [11,22,23]. This fact emphasizes that for diffusion over closed surfaces, the observed motility depends on the curvature of the surface and on the plane of observation [24][25][26].
We present the results of the MSD for interacting particles in a crowded spherical surface by using Molecular Dynamics (MD) simulations and implementing five different pairwise-interaction potentials between the particles in a relatively high surface density medium (ρ p 0.92). The particles are free to move on the entire spherical surface of radius R. To keep the particles on the surface, we use the algorithm described by Juffer et al. [27]. The continuous interaction potentials used are different cases of the generalized form of the Mie potentials [28]: LJTS, WCA(12.6), WCA (50.49). Besides, some simple discontinuous potentials like the shoulder-square (SS) and the well-square (SW) were adapted to their continuous-mathematical form described by Munguía-Valadez et al. 1 With the aim to provide a physical interpretation of the results of these numerical calculations, we formulate an analytic model based on the generalized Smoluchowski equation with timedependent coefficients [29][30][31][32][33][34]. Using this tool, we can identify how the change of perspective, confinement and crowding are incorporated in an effective diffusion coefficient that takes into account the entropic confinement and the anomalous diffusion. The model used here was derived in the context of mesoscopic non-equilibrium thermodynamics (MNET) and has been widely used in order to obtain kinetic equations for transport phenomena, like diffusion-adsorption processes, anomalous diffusion, activated processes, diffusion in pores, and diffusion in the presence of entropic barriers [35][36][37][38][39]. In particular for our purposes, MNET has been successful in describing diffusion on other confined systems [38,[40][41][42].
The comparison between our numerical simulations and thermodynamic-based model allows us to propose a new interpretation of the observed dynamics of the MSD for diffusion of interacting particles over a sphere, and to compare the different subdiffusive regimes associated with the hardness or softness of the particle collisions. We have chosen a spherical geometry in order to rule out the effects of different local curvatures in the surface. This will also allow us to recover previous theoretical results deduced for the sphere in absence of interaction.
Regarding the organization of work, in Section 2 we present our Molecular Dynamics simulations on the diffusion of finite sized-particles over an spherical surface for different interaction potentials among particles. Then, we study the dynamics of the diffusion and the behavior of the MSD using the results provided by a Smoluchowski description which is presented in Sections 3.1 and 3.2, for free and interacting particles, respectively. Comparison between simulations and model is presented in Section 4. Finally, discussion and concluding remarks are provided in Section 5.

MOLECULAR DYNAMICS SIMULATIONS
We present the details of the simulations as well as a brief description of the interaction potentials used in this work.

Interaction Potentials
The continuous interaction potentials are defined through: where ϕ n,m (r ij ) is the most general form of the Mie potential [28]. For example, we can write the Truncated and Shifted Lennard-Jones (LJTS) interaction by choosing n, m 12, 6, and then ϕ LJTS r ij ϕ 12,6 r ij − ϕ 12,6 (r c ), r ij ≤ r c , 0, where σ( 2a) is the particle diameter, r ij r i − r j the distance between the centers of mass of the i-th and the j-th particles, ϵ the potential well depth, and r c ( 2.5σ) the cutoff radius. On the other hand, a generalization of the Weeks-Chandler-Anderson potential (WCA) [43,44] can be written as where r ij (n/m) 1/n− m σ is the value that corresponds to the minimum potential. In our case, we use the parameters n, m 12, 6, [WCA (12,6)] and n, m 50, 49, [WCA (50,49)] to incorporate pseudo hard-sphere interaction models into the study, see Figure 1B [44]. The Square-Well (SW) and Square-Shoulder (SS) potentials can be approximated with the Generalized Continuous Multiple Step Potential (GCMS), which incorporates a contribution of excluded volume followed by multiple steps that model repulsive barriers or attractive wells as the case may be. Its simplified form can be written as: where ω defines the spatial extent of the particle core in units of σ, q is related to the hardness of the potential, and a 0 is a factor that defines an attractive step (a 0 < 0) or a repulsive one (a 0 > 0). The parameters used to reproduce the SS and SW are: ω 0.5, q 500, a 0 1 and a 0 −1, respectively, see Figure 1C y Figure 1D.

Simulation Method
In this work all quantities are assumed to be expressed in conventional reduced units, with m, σ and ϵ as the units of mass (set equal to 1), distance, and energy, respectively. According to this convention, the temperature (T p ) is in units of k B T/ϵ where k B is Boltzmanns constant, the density (ρ p ) in units of ρσ 2 , the time (δt p ) in units of δt(ϵ/mσ 2 ) 1/2 , and the energy (ϕ p ) in units of ϕ/ϵ. A monodisperse system of spherical particles was studied with the restriction of moving on the surface of a sphere of radius R 9.3σ. To study the dynamic behavior of each system, Molecular Dynamics simulations were performed in the canonical assembly (NVT) for a total of N 1000 particles embedded in the spherical surface, with a surface density given by ρ p ( N/A) 0.92. The R and ρ p values were taken for compatibility with previous studies by Vest et al. [45][46][47]. The system was placed in a thermal bath at a temperature T p 1.0 set constant by using the isokinetic method; the integration of the equations of motion was carried out using the velocity-Verlet algorithm [48], with the restriction f i (r i ) |r i | 2 − R 2 0 (see Ref. [27]) and a time-step δt p 10 − 3 for LJTS, WCA (50,49) and WCA(12, 6) systems, and with δt p 10 − 4 for SS and SW systems. The simulations were run by 10 6 time-steps for LJTS, WCA (50,49) and WCA (12,6), and by 10 7 time-steps for SS and SW to reach thermodynamic equilibrium, then for 3 × 10 7 and 4 × 10 8 time-steps, respectively, for the calculation of the mean square displacement. Some representative trajectories of a tracer particle for each potential are plotted in Figure 2.
The results of our simulations for the MSD are plotted in Figure 3. At the right side we plot the geodesic MSD over the sphere with the saturation valor of (π 2 − 4)R 2 /2. At the left, we plot the MSD measured when the displacements of the particles are projected the plane of observation XY, giving the saturation value of 2R 2 /3. These saturation values for long times only depend on the geometric configuration and therefore are the same for free and interacting diffusing particles [10].
The MSD plotted in Figure 3 shows four different stages: 1) the ballistic regime where the MSD increases as t 2 , 2) a transition between ballistic and diffusive regime where the slope of the MSD decreases and changes of concavity, 3) the subdiffusive regime where, as we will calculate below, the MSD increases as t α with α(1 and, finally, for long times, 4) the saturation regime where the MSD no longer increases since the domain of the sphere is finite. These four regimes are seen in the geodesic and projected version of the MSD plotted at right and left of Figure 3, respectively. Some or all of these four stages have been found both in numerical simulations [9,[49][50][51][52] and experiments, mainly for protein and lipid diffusion on cell membranes [53][54][55][56].
The form of the plots for the MSD in Figure 3 depends on the combined effect of the curvature of the surface, the interaction of the particles and, in the case of the projected displacement, the projection on the observation plane. In the next section we will deal with the connection between these different effects and their effect in the measured effective diffusion coefficient.

THE MSD OF FREE AND INTERACTING PARTICLES
Commonly, experimental setups for measuring the diffusion of particles on curved surfaces use confocal microscopy techniques [33,57,58]. Due to the focal length inherent to these techniques, the measure of the motions of the particles under study is preferentially performed in terms of their projection to the focal plane [22,23]. In view of this, in the present section we study the equivalence between the description of the free diffusion over the surface (the geodesic displacement) and compare with the evolution of its projection in the planar disk constrained by the sphere, see Figure 4. First, we will deduce an expression for the MSD of free particles moving at a flat disk by considering how the entropic restriction of the movement (imposed by the fact that the particles cannot travel beyond the radius of the disk) is reflected in the transport properties that are used to describe the motion. Then, we will modify the MNET description to include the interaction and crowding of particles.

The MSD of Free Punctual Particles
From a thermodynamic point of view, the problem can be tackled by estimating the equivalent force keeping the system confined to diffuse in the planar disk delimited by the surface of the sphere of radius R, see Figure 4. This equivalent force can be written in terms of a Taylor expansion. To first order, this correspond to an harmonic-like potential of the form The potential here is of an entropic origin, since it is due to the fact that the particles cannot leave the sphere projection or, equivalently, the planar disk. The maximum value of this potential, associated with the equivalent force, can be estimated by using the equipartition theorem [59]. Comparing the mean kinetic energy of the particles, K (3/2)kT (which corresponds to one half for each degree of freedom in the three dimensional description), with the maximum potential energy available for the trapped particles, U max (1/2)κ 0 R 2 , we can estimate the value of the effective restorative coefficient κ 0 as Hence, the projected motion of the particles constrained to diffuse in a disk under the influence of the equivalent harmonic force can be described by means of the Smoluchowsky equation where k B and T are Boltzmann constant and temperature, respectively. Since the equivalent force has radial symmetry, assuming an initial configuration with the same symmetry allow us to reduce the calculation of the MSD to the contribution in the radial coordinate r. In this case, we have where the force that constraint the particles motions is F −∇U, with U given by Eq. 5, and the gradient and Laplacian operators are represented in polar coordinates.
The mean square displacement in the constrained disk system is therefore given by The time evolution of the MSD is obtained by taking the derivative And substituting Eq. 8. Assuming zero flux boundary condition: zp/zr(R) 0 in the border of the disk and the conservation normalization equation, R 0 prdr 1, we can prove, using successive integration by parts, that The comparison of the previous result with the result of diffusion in an infinite planar surface, 〈r 2 〉 4D 0 t, shows that the effect of the confinement is proportional to the ratio of the potential and kinetic energies: where the average here is over the ensemble of non-interacting particles. Eq. 11 is readily solved using that 〈r 2 〉 → 0 when t → 0 as Substituting the estimated value of κ 0 given in Eq. 6, valid when particle interactions are negligible, we obtain This result has been obtained by using different methods in Ref. [10] and, as expected, only depends upon geometrical considerations resulting from averaging methodologies applied directly to the Laplace-Beltrami diffusion equation.
This result describing the projection of the surface to the disk can be contrasted to the result of particles really moving on a planar disk of radius R, where K k B T (since there are only two freedom degrees), κ 0 2k B T/R 2 , and therefore, from Eq. 11, we have Comparison between Eqs 14, 15 shows that the displacement in a sphere respect to a planar projection has an average reduction by a factor of 2/3 for long times. This reduction factor in the displacements can be deduced also by evaluating the ratio between the effective diffusion coefficient measured on a planar region, D 0 , with respect to that measured in the observation plane, D ⊥ . The result is Where 〈n 2 z 〉 is the geometric average of the quadratic vertical component of the normal vector [10,59]. For the sphere we have n z cos θ, from where it is straightforward deduced that 〈n 2 z 〉 1/3 giving the expected result D ⊥ /D 0 2/3.

The MSD of Interacting Non-punctual Particles
In the previous subsection we have shown that the entropic restrictions, present by the fact that tracer particles move over a surface with a finite number of accessible states, are responsible for the observed saturation value reached by the MSD at long times. However, as we have seen in Figure 3, the simulations show that the time behavior of the MSD at intermediate times scales with a time dependence t α with α ≤ 1. Therefore, for certain interaction potentials, anomalous diffusion is observed. This fact requires that the model developed in Section 3.1 should be generalized to cope with particle interaction.
This generalization goes along the same lines of Santamaría-Holek et al. [32], where a Smoluchowski description with FIGURE 4 | Projection on an equatorial plane of the Brownian motion which takes place over the surface of a sphere of radius R. The polar radial coordinate of the particle r is measured on the plane.
Frontiers in Physics | www.frontiersin.org March 2021 | Volume 9 | Article 634792 time-dependent coefficients was used to demonstrate that the exponent characterizing the sub-diffusion is controlled by the nature of the local cages and the free space at disposal for their motion. In our case, since the single particle MSD reported here accounts for the motion of independent particles in a crowded media, the description of the diffusion can be given also in terms of the generalized Smoluchowski equation for the single-particle distribution function [61]. This equation for the probability is the counterpart of the respective generalized Langevin equation with memory effects, commonly used to describe microrheological experiments and has the peculiarity of being equivalent in the long-time regime to the Smoluchowski equation with time dependent coefficients [29,30]. This fact was generalized in Refs. [31,33,62] for the case of power-law dependent memory kernels which are equivalent to time dependent diffusion coefficients of the form where τ is a characteristic time of the anomalous diffusion process. In this equation, the exponent α ∈ (0, 1] characterizes how strong sub-diffusion is and, as it was shown in Ref. [32], has the general form Besides the entropic effect introduced through the ratio k B T/κ 0 R 2 , the Eq. 18, through the factor 1 − B 1 (a/ξ), contains the two main ingredients mentioned previously, namely, the freearea at disposal to perform diffusion a/ξ, and the nature of the local cages through the coefficient B 1 . The characteristic length ξ is related with the free-area in which the finite particles can move in a surface of total area A: ξ A 1/2 free xR(1 − Ma 2 /R 2 ) 1/2 where M is the number of particles forming the cage. In this case, we have This expression for the sub-diffusion exponent shows that it only depends in the radii of the particles and sphere, and on the parameter B 1 . This parameter indirectly depends in the nature of the interaction potentials used. In the context of hydrodynamics [63], it was related with the correction introduced by hydrodynamic interactions over the motion of a particle when this motion takes places near to a solid wall. Since these interactions are related to the potential interactions among particles, one may conjecture that the parameter B 1 is a measure of the effect that different interaction potentials have on the cages' dynamic structure and, therefore, of the different values that the sub-diffusion exponent may take.
Even in this approximation, it results difficult to estimate theoretically the parameters in Eq. 18. However, it is clear that the exponent of the subdiffusive process is linked to the shape of the interacting potential, the curvature of the surface and the projection of the forces to the observation plane. In order to find its value for the different potentials in Figure 1, we will adjust the data of the numerical simulations, using the Smoluchowski equation, but now in terms of this time dependent diffusion coefficient as we detail below.
Therefore, we will consider that the single particle Smoluchowski equation describing the dynamics of a tracer particle in an effective medium of N − 1 interacting nonpunctual particles takes the form zp zt Where now the time dependent coefficient is given through Eqs. 17,19. From the last equation, repeating the same procedure of the last section, the temporal evolution of the MSD is obtained to be Using Eq. 6 we finally have Which is the general expression for the MSD of the tracer particles valid for the whole time interval. From this equation, it is possible to deduce the temporal behavior of the diffusion coefficient for all the stages of the process in terms of the numerical data provided by the simulations as .
Notice that this equation reproduces the expected behavior for short and long times. For short times 〈r 2 〉 → 0 and d〈r 2 〉/dt → 4D 0 and therefore D s → D 0 . For long times, it is expected that the MSD saturates and, therefore, d〈r 2 〉/dt → 0 with D s → 0.
At the other hand, when (Eq. 22) is used in combination with the diffusion coefficient in Eq. 17, it accounts only for the process once the sub-diffusion stage has started. At intermediary times, it accounts for the sub-diffusion associated with particle interactions and cage diffusion effects and, for long times, it account for the fact that the particles have a restricted number of accessible states in a finite domain. In the next section we will use the numerical simulation in order to find the general behavior of the temporal diffusion coefficient and the exponents of subdiffusion for the potentials we have considered in Section 2.

COMPARISON BETWEEN MODEL AND SIMULATIONS
We have shown that the projection on the equatorial plane of the geodesic MSD, directly measured in the numerical simulations, is an equivalent form to characterize the diffusive dynamics of the particles. We have explained how the saturation value of the MSD rescales due to the restricted number of accessible states, giving rise to entropic constrictions on the dynamics of the particles.
For short and intermediate times, our methodology allows us to consider how different aspects modify the behavior of the MSD of interacting particles forced to diffuse in the surface of the sphere. All these aspects are clearly represented in Figure 5, where we plot the time-dependent diffusion coefficient, D s , given in Eq. 23 for the potentials chosen, and illustrated in Figure 1. For constructing this plot, the data from the functions 〈r 2 〉(t) and d〈r 2 〉/dt obtained from numerical simulations are directly substituted in (Eq. 17) for D s (t).
The first aspect to notice is the reference diffusion coefficient, D 0 . This is marked as a dashed line in Figure 5 and represents the surface diffusion coefficient in an infinite planar surface, and negligible interaction potential among particles. In Figure 5 it is seen that, for the used potentials, the value of the diffusion coefficients augments according to the order: SW, LJTS, SS, and WCA (12,6). The color key is the same as in Figure 1.
The second aspect to notice in the temporal dependence of the diffusion coefficient D s (t) is the emergence of sub-diffusion due to the crowding effects. This results in different slopes of the curve D s (t) at intermediary times and, therefore, corresponds to different values of the sub-diffusion exponent α. In order to compare quantitatively the sub-diffusion for the different potentials, in Figure 6, we adjust the data from numerical simulations using (Eq. 22) with the temporal dependence of the diffusion coefficient in (Eq. 17) using the least squares method. The value of the local diffusion coefficient D 0 and the sub-diffusion exponent α together with the fits of the data are given in Figure 6. Figure 5, show that the value of the diffusion constant D 0 is smaller for hardest potentials (WCA (12,6) and SS) than for softer potentials (WCA (50,49), LJTS and SW). The potentials having wells also show this tendency since the diffusion constant for the LJTS is larger that of the SW potential. The attractive part of the potentials decreases, in turn, the kinetic energy of the particles surrounding the tracer particle and, therefore, the energy availability to perform position fluctuations. This interpretation is clear after a comparison between the potentials WCA (50,49) and WCA (12,6), which have a similar functional form but different degree of penetrability. Another aspect to be noticed is the steepest drop of the diffusion coefficients (lower α) for the potentials with a tilted repulsive part after r 1, that is, SS and WCA (12,60). A better understanding of the effect of the shape and hardness of the interaction potentials on the diffusion constants and the subdiffusion process will require its own work.
Finally, the third aspect to notice is that both the ballistic and long-time behaviors of the MSD appear in all cases since not depend on the interaction potential. Because of this, they do not provide information of the fluid medium where the diffusion occurs. As it is seen in Figure 5, the final drop of the diffusion coefficient to zero occurs earlier as the diffusion coefficient is lower, as expected.
Our work shows that the deviations of the MSD relative to the planar behavior represent the effect that the particle interaction and crowding effects have on the dynamics. This effect is coupled with the presence of an entropic force restricting the number of accessible states. We prove that these two effects are well captured by the projected diffusion coefficient D s , given in Figure 5 and which has been reported in previous numerical simulations and experiments [53,64].
The excellent agreement in Figure 6 among simulation results and the theory proposed indicates that particle interactions control the cage effect and the sub-diffusion regime. The delay induced by the crowded scenario appears since the transient stage, and its influence is captured by the memory effects that the effective diffusion coefficient D s (t) incorporates. We have found that the value of the parameters used for fitting the MSD curves depend upon the used potentials and further numerical and theoretical studies are necessary in order to show the relation between the form of the potential interaction and the surface diffusion coefficient. Similar commentaries can be done concerning the specific value for the radii of the particles and the surface density. These values allowed us to capture the effect of the surface density and the form of the potential on the subdiffusion observed. Nonetheless, a more detailed study on the dependence of the diffusion properties on this factor will require a future work.

CONCLUSION
We have studied the statistical and dynamical properties of a system of interacting particles constituting a fluid embedded on a sphere's surface. To understand the link between the microscopic properties of the fluid and the observed MSD, we have proposed a diffusion model incorporating memory and entropic effects that allowed us to derive a simple and powerful expression for the projected time-dependent diffusion coefficient. This expression allows us to characterize the dependence of the diffusive dynamics on the particle interactions and surface curvature. This approach takes into account the change of perspective inherent to the experimental measure of diffusion of particles on closed surfaces. The formalism proposed in this work reproduces also some limiting cases previously studied in the literature and exemplifies, with the use of ideas emerging from irreversible thermodynamics, how we can provide a more flexible description that helps to interpret complex phenomena occurring at the surface.
In this work, we have shown in Figure 3 that the dynamics of the MSD reflects four stages. Two, at short times t p ≤ 10 − 1 s, are associated with the ballistic regime and how it ceases; and other two (t p ≥ 10 − 1 s) associated with the diffusion regime in presence of cages at the intermediate and long time regimes, corresponding to anomalous diffusion and saturation, respectively. Using the ideas of the MNET we use the Smoluchowski equation for describing the process once the diffusion regime starts and this is shown Figure 6 for the intermediary and long time behaviors of the MSD. In contrast, modeling the regimes at short times will require another kind of kinetic mechanism. However, it is worth mentioning that, under the perspective of biological applications of surface diffusion in crowded membrane environments we discuss below, only the subdiffusive regime is accessible to measurements with the actual experimental timeresolution capabilities and therefore, this does not constitute a limitation to our conclusions.
In the realm of applications, this study allows us to establish a connection with models and experiments studying the diffusion of several tracers in the cytoplasm of living cells which often exhibit heterogeneous distribution of macromolecular crowders [65][66][67]. This crowding is extremely important in the surface membrane and might affect the properties of anomalous diffusion. It is known, for example, that going from a less to a more crowded region will slow down the dynamics, both in terms of the exponent and diffusivity [68][69][70]. This results very important since the slow transport in the cell membrane of lipids and proteins is linked with protein cluster formation, phase segregation, lipid droplet formation, signal propagation and other crucial functions occurring on the cell surface which can be enhanced by the presence of sub-diffusion [71,72].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.