An efficient approach to study membrane nano-inclusions: from the complex biological world to a simple representation

Membrane nano-inclusions (NIs) are of great interest in biophysics, materials science, nanotechnology, and medicine. We hypothesized that the NIs within a biological membrane bilayer interact via a simple and efficient interaction potential, inspired by previous experimental and theoretical work. This interaction implicitly treats the membrane lipids but takes into account its effect on the NIs micro-arrangement. Thus, the study of the NIs is simplified to a two-dimensional colloidal system with implicit solvent. We calculated the structural properties from Molecular Dynamics simulations (MD), and we developed a Scaling Theory to discuss their behavior. We determined the thermal properties through potential energy per NI and pressure, and we discussed their variation as a function of the NIs number density. We performed a detailed study of the NIs dynamics using two approaches, MD simulations, and Dynamics Theory. We identified two characteristic values of number density, namely a critical number density nc = 3.67 × 10−3 Å−2 corresponded to the apparition of chain-like structures along with the liquid dispersed structure and the gelation number density ng = 8.40 × 10−3 Å−2 corresponded to the jamming state. We showed that the aggregation structure of NIs is of fractal dimension dF < 2. Also, we identified three diffusion regimes of membrane NIs, namely, normal for n < nc, subdiffusive for nc ≤ n < ng, and blocked for n ≥ ng. Thus, this paper proposes a simple and effective approach for studying the physical properties of membrane NIs. In particular, our results identify scaling exponents related to the microstructure and dynamics of membrane NIs.


Introduction
A lipid bilayer is a thin polar membrane composed of two layers of lipid molecules. These membranes form a continuous barrier around cells and are essential for their homeostasis by regulating the ions' and molecules' diffusion through them. 1,2 The cell membrane of almost all living organisms and many viruses are made of a lipid bilayer, 3 just like the membranes surrounding the cell nucleus and organelles. The lipid bilayers are impermeable to hydrophilic molecules and particularly to ions. This selectivity enables the cells to regulate, in particular, the pH and salinity of their cytosol by using transmembrane proteins. Thus, provide a transporter function capable of generating and maintaining a concentration gradient and the extracellular environment. 4 Biological membranes usually contain other constitutions than phospholipids. These proteins and other integrated components are involved in several pathologies, including cancer. 5 In this context, W. Yang et al. reported a mechanism by which the anti-tumor response of cells used in cancer immunotherapy can be potentiated by modulating cholesterol metabolism. 6 The cholesterol level in cell membranes determines the physicochemical properties of these membranes. 7 In addition to these organic constituents, the membrane may contain inorganic nano-inclusions (NIs), such as small and macro-ions, or more complex structures that can play important roles. In this context, A. O. Elzoghby et al. show that protein-inorganic hybrid nano-inclusions offer promising opportunities for drug delivery and tissue imaging, particularly in the eld of oncology. 8 Besides, they are the lifeblood in thin-lm composite membranes used in desalination. 9 NIs impose a well-dened membrane deformation is the subject of intensive studies due to the potential opportunities that they offer for biomedical applications in diagnosis and therapy because of their rigidity and their perfect monodispersed character. [10][11][12][13][14][15] To our knowledge, the techniques for determining the physical properties of nano-inclusions remain rare and still to this day a promising eld to explore, in particular using MD simulation. However, many theoretical and experimental studies have been interested in determining the interactions between NIs across the membrane. [15][16][17][18][19][20] In this context, I. Koltover and coworkers presented experimental studies of the aggregation of spherical colloidal particles adhering to lipid membranes. 18 They observed that two-particles aggregation was found consistent with a short-range-attractive interaction and no obtained sign of a long-range one. Surprisingly, triplets of particles form chains, and a ring-shaped aggregate is observed around the vesicle size. Also, linear chain arrangements are obtained in simulations of similar situations for some particle size and adhesion regimes, 19,20 using a scaling argument to show that this was not due to membrane-mediated interactions. 20 But it is due to the adhesion of the particles to the membrane since a linear aggregate gives a higher adhesion area than a compact one. Therefore, such a phenomenon would not occur in the case of inclusions. 16 Costanza Montis et al. performed several experiments to compare the interaction of Extracellular Vesicles-derived (EVSLB) and synthetic Supported Lipid Bilayers (SLBs) with cationic superparamagnetic iron oxide nanoparticles (SPIONs). Their experimental results revealed that the SPIONs-EVSLBs interaction is signicant in comparison to SPIONs-SLBs. 17 Particles that form jammed structures are ubiquitous in nature and omnipresent in many consumer products. 21 Being metastable solids, they are mechanically rigid structures that generally form at low particle density due to a spatial cluster of aggregated particles. The colloid aggregation is controlled by the number density, the attractive interaction potential, and the temperature. 22 Colloid aggregations in spanning network systems are known as gelation or coagulation. 23,24 Thus, several works have been devoted to determinating the physical parameters and demonstrating phase boundary mapping of the phase transition by examining a range of volume fractions, temperature, and interaction strength between particles. [25][26][27][28][29][30][31][32][33][34][35] As schematically shown at the top of Fig. 1, membrane nano inclusions interact with the surrounding membrane lipids. Also, they interact with other inclusions. Furthermore, they can form complexes with other membranes NIs and with integral proteins. Thus, an all-atom MD simulation has many details that cost a lot of time in computing and programming skills. A simple alternative is to determine the direct interaction between the particles under study, where the other particles are eliminated and replaced by their effects. These effects appear in the expression of the simplied interaction potential with a relatively simple mathematical formula. From a mesoscopic viewpoint, we consider the membrane constituents as a twodimensional homogeneous uid since these biological constituents have almost the same structural and dynamics properties. 36 In this case, we simplify the study of the structure and dynamics of NIs into a two-dimensional colloidal system with an implicit solvent (membrane). D. Constantin and coworkers, in 2008, determined the direct interaction potential between BuSn12 NIs inserted within diluted surfactant bilayers by using X-ray synchrotron scattering. 15,37 Nine years later, M. Benhamou and coworkers added a prominent correction to describe the interaction between NIs within biological membrane bilayer. 38 In the work of Constantin et al., 15,37 the repulsive interaction between inorganic BuSn12 nanoparticles inserted into surfactant bilayers was measured using Small Xray Synchrotron Scattering. Accordingly, the analysis of the measured structure factor of the two-dimensional uid formed by the inserted nanoparticles led to their membrane-mediated interaction. This interaction is repulsive and well tted to a Gaussian function. Thus, M. Benhamou et al. 38 have added an attractive term to the repulsive interaction form proposed by D.
Constantin. This form can describe will the NIs interaction within a biological membrane bilayer.
In this case, we simplify the study of the structure and dynamics of NIs into a two-dimensional colloidal system with an implicit solvent (membrane). From previous works, D. Constantin and coworkers, in 2008, determined the direct interaction potential between BuSn12 NIs inserted within diluted surfactant bilayers by using X-ray synchrotron scattering. 15,37 Nine years later, M. Benhamou and coworkers added a prominent correction to describe the interaction between NIs within biological membrane bilayer. 38 In the work of Constantin et al., 15,37 the repulsive interaction between inorganic BuSn12 nanoparticles inserted into surfactant bilayers was measured using Small X-ray Synchrotron Scattering. Accordingly, the analysis of the measured structure factor of the twodimensional uid formed by the inserted nanoparticles led to their membrane-mediated interaction. This interaction is repulsive and well tted to a Gaussian function. Thus, M. Benhamou et al. 38 have added an attractive term to the repulsive interaction form proposed by D. Constantin. This form can describe will the NIs interaction within a biological membrane bilayer.
In this work, we expect the NIs to undergo a phase transition, corresponding to the critical and the gelation numberdensities, that we note n c and n g , respectively. For these densities, we expect the physical properties to present In the MD simulations, we considered the membrane as a two-dimensional implicit liquid for the sake of simplicity. Thus the nano-inclusions are treated as a two-dimensional colloidal suspension interacting via an efficient interaction U eff . Then, these considerations reduce several details and then allow the study of a large number of nanoinclusions. The numerical data of the radial-distribution-function g(r), the mean-square-displacement MSD, and the diffusion coefficient D c allow the determination of the gelation number density n g , the fractal dimension d F , and other related parameters. a discontinuity. For examining this hypothesis, we apply MD simulations, Scaling Theory, and Dynamics Theory to investigate the structure and dynamics of NIs in bilayer membranes. For this, we shall perform a series of simulations for various number densities. Thus, we determine n c and n g , and we quantify the relevant scaling exponents. For this purpose, the structure and dynamics behavior will be discussed using a detailed analysis of the radial distribution function and the mean square displacement. This paper proposes a simple and effective approach for detecting the gelation phenomenon, relevant for many aggregation processes at the level of biological membranes. In particular, our results identify scaling exponents related to the microstructure and dynamics of membrane NIs.

Interaction potential
In this work, we considered an assembly of monodisperse spherical inorganic nanoparticles immersed in a twodimensional bilayer-membrane. We assumed that these nanoinclusions are rigid and charge-neutral. We denoted by n ¼ N S and d ¼ 2R, the nanoparticle number density, and their common diameter, respectively. N accounts for their number, S for the membrane area, and R stands for their radius. Because of the membrane undulations, the immersed nanoparticles experience a membrane-mediated repulsive force. The latter was measured using Synchrotron Small-Angle X-ray Scattering. 15,37 The interaction potential determined with the help of the structure factor measured using RPA. 39 The latter compared to the theoretical one, which explicitly contains the pairpotential. The reference was the structure factor of a hardcore disc model. Herein, we recall simply the expression of the membrane-mediated repulsive potential, where U R (r) represents the repulsive part of the total interaction potential, and r is the separation distance between the particles. The potential amplitude A, whose magnitude is in the order of some thermal energies k B T. r c is a cutoff distance, for which the interaction potential tends to zero, used in the MD simulations to reduce the calculation time, and z represents the potential range. We give an explication to this distance-scale by saying that it presents the length of the correlation in the plane (size of the hole and the valley of membrane's undulations). The latter depends on the bilayer membrane characteristics 40-42 through the modulus of curvature, k, and the interfacial tension, g. According to M. Benhamou et al. corrections, 38 the characteristics of the membrane k, g, the potential amplitude A(n) and the renormalized range z(n) depend on the number density of NIs, n. Thus, to propose an explicit variation-law to present these quantities as a function of the density, the authors assert it is a typical value, n 0 , where the bilayer membrane ceases to be uid, then becomes rigid. Therefore, at n ¼ n 0 , the ripple of the membrane disappears, that is, z(n) ¼ 0. This implies that A(n) ¼ 0. Also, the authors say that the special value n 0 corresponds to the situation where the nanoparticles are extremely congested, then cover the entire surface of the membrane. Thus, n 0 is dened as follows: n 0 ¼ 1 pR 2 , where pR 2 is the disc surface. Thus, they assume, intuitively, that the amplitude A(n) and the range z(n) must decrease when the density n increases. For these quantities, they expect the following approximate forms, here, z 0 denotes the bulk value of the in-plane correlation length (in the absence of nanoparticles). The constant A 0 that measures the number of contacts between interacting NIs should be proportional to the inverse of the squared bending modulus k, 43-45 that is A 0 $ k À2 , where the proportionality coefficient is known [43][44][45] and presents as a combination of the thermal energy k B T and the nanoparticle diameter, d. Therefore, as it should be, the latter vanishes for the rigid bilayermembranes k / N. Besides the above repulsive potential, the NIs experience naturally an attractive force of van der Waals type, 46 whose corresponding potential reads (in Derjaguin approximation 47 ), where d ¼ 2R is the NIs diameter and A H denotes the Hamaker constant, which is in the 10 À22 to 10 À19 J range. Therefore, the overall effective potential U eff (r) is, Thus, the pair interaction potential between membrane NIs is the sum of a repulsive term, originating from the bilayermembrane undulations, and an attractive one derived from the polarity of the host medium. We note that, in contrast to the repulsion, the VDW attraction is negligible for low-density NIs, but it is indispensable for high-density NIs.

Simulation details
In this work, we discuss a particular type of spherical NIs, namely the BuSn12 nanoparticles. Where BuSn12 denotes {(BuSn) 12 O 14 (OH) 6 } 2+ (4-CH 3 C 6 H 4 SO 3 À ) 2 , the butyltin oxo cluster. The synthesis and characterization of BuSn12 is described in detail by Eychenne-Baron et al. 48 BuSn12 inclusions are considered as a spherical nanoparticles of radius R ¼ 4.5Å and molecular weight of 2866.7. The corresponding parameters of the effective interaction potential, discussed above, are as follows: In the come aer, and for a purpose to optimize the calculation time, we use the dimensionless unit, meaning that all quantities are unitless in the MD simulations. 49,50 Accordingly, the distance is reduced by the NIs radius R, the number density is reduced by R À2 the energy is reduced by k B T, the temperature is reduced by the room temperature T a ¼ 298.15 K, the time is reduced by , and the pressure is reduced by k B T R 2 . In this work, we calculated all dynamics and structural properties using NST MD simulations by analyzing the trajectory of each particle. The simulated system consists of 10 6 colloidal particles in a square lattice, with applied periodic boundary conditions in the two dimensions. The NST statistical ensemble conditions are adjusted using temporal integration on non-Hamiltonian Nose-Hoover motion equations. These equations generate the positions and the velocities sampled from the canonical ensemble (NST). Thus, adapt these kinetic quantities of NIs at each time step. This style xes the temperature by adding dynamics variables coupled with particle velocities. The equations of motion used are those of Shinoda et al., 51 which associate the hydrostatic equations of Martyna et al. 52 with the strain energy proposed by Parrinello and Rahman. 53 Temporal integration schemes closely follow the Verlet integrators, preserving the measures and reversible in time, derived from Tuckerman et al. 54 We performed the MD simulations using the LAMMPS soware package. 49 We started from a random initial conguration. Then, the initial conguration is equilibrated in the NVT ensemble before the production of the results. In this way, any proposed initial arrangement leads to the same equilibrium. We run the simulations at different densities above and beyond the gelation density n * g . The system is balanced long enough (1 000 000 time-step) with a time-step Dt* ¼ 0.01.
To analytically determine n c and n g , we go back to the effective interaction potential. We note that, in addition to the number density, there are other relevant parameters that we have considered xed, namely the temperature and the NIs size. Practically, we characterized U eff by two behaviours, attractive or repulsive, in other words, negative or positive. From the literature, M. Benhamou et al. found that n c is unique and given by: 38 with a a dimensionless coefficient The critical number density, n c , is of the order of n c ¼ 3.6 Â 10 À3 A À2 , for BuSn12 NIs with the effective interaction potential parameters: From the theoretical investigation, we summarize the following number density regimes comparing the number density n with the predicted threshold of n c : (1) Low-density regime corresponding to n < n c : in this regime, the potential exhibits two zeros and two extremums, with an energy barrier between the two zeros. The rst extremum is a minimum around which the immersed nanoparticles occulate. Therefore, in the low-density regime, the potential presents a high energy barrier that prevents NIs gelation.
(2) First critical density regime corresponding to n z n c : for this critical value, the potential is always attractive and exhibits one maximum, which is a zero at the same time for r ¼ z(n c ) and a secondary minimum at the right. At this specic value of density, the energy barrier disappears.
(3) High-density regime corresponding to n > n c : in this regime, the potential is purely attractive and presents as an increasing monotone function of the surface-distance r. This is the region of a prevised gelation number density n g . The rst gelation number density will be determined using MD simulation.

Structural properties
The structural analysis has been conducted using the twodimensional Radial Distribution Function (RDF) and Coordination Number (N c ).
The RDF is dened as the probability of nding two NIs at a distance r, each to other.
where N is the number of NIs and S is the area of the host membrane. In the MD framework, we performed the above sum over all pairs of nano-inclusions. The average density distribution of NIs, at any point, is referred to as the bulk-density, n b . Thus, the average density of the NIs at a given distance, r, from the center of any considered NIs, n r , is related to RDF by: We dene the coordination number N c as the number of neighbors NIs that are within the specied cutoff distance of the pair interaction potential from a central inclusion, where r c stands for the position in which the interaction potential between NIs goes to zero. Thus, the running coordination number, N(r), is dened by, herein, r is an arbitrary distance from the center of the tagged nano-inclusion.
3.1.1 Scaling theory predictions. As the density increases, the attractive forces begin to dominate 23,24 and lead to a macroscopic aggregate of NIs. As a consequence, inclusionsdensity increases locally. In practice, a simple descriptive denition of the fractal dimension d F is adopted, 55-58 we consider a disc of radius r centred at some nano-inclusion taken as the origin of space, and we estimate the number of adjacent NIs inside it. Since the distribution of the NIs is random, the disc (or the aggregation space) is of fractal structure, with a fractal dimension denoted as d F . The latter is naturally less than the membrane surface dimension (d ¼ 2). As a result, the number of NIs inside the disc that is directly proportional to their total mass identies with the current coordination number, N(r), and we write the scaling relation, Using the denition of N(r) and its relation with RDF mentioned in the above section (eqn (12)) suggests that RDF behaves, as a function of the separation distance, as: with H a positive dimensionless amplitude, s ¼ 2R accounts for the NIs diameter. Since d F < 2, RDF decreases with distance, by the fact that the correlations between the NIs attenuate.
The rst precious information we extracted from the RDF scaling-law deals with the principal peak height g max and its width r max variation as a function of the NIs number density. Such a law is also valid for r ¼ r max , and we then write, or in log-log scale, Then, the slope of the curve representing ln(g(r max )) as a function of ln(r max ) makes it possible to determine the fractal dimension of the NIs aggregate, d F . The Scaling Theory suggests that d F does not depend on the number density n. Thus, the Scaling Theory predicts the d F to be a universal quantity. The amplitude, H, can be obtained by extrapolation on the vertical axis. We recall that since, as part of the Scaling Theory, d F < 2, g max decreases with the position of the rst peak of RDF, r max .
3.1.2 Structure from MD simulation. Fig. 2 depicts instantaneous snapshots at different number densities at T ¼ 298.15 K. For low number density, more free-area is available, and the correlation between the NIs is weak. Analysis of the equilibrated particles conguration showed the apparition of small chains-like structures. Thus, these form small barriers and slightly limit the NIs diffusivity. At the critical number density n c , NIs form large aggregates manifesting as percolated gels (chains-like) and jammed clusters. Thus, each cluster gets stuck in a cage formed by its neighbouring clusters. In this case, the correlation between NIs is signicant. Thereby, we expect slow dynamics. At the gelation number density n g , the NIs form crowded chains-like structures. In this case, the particles cannot move freely because of the dominant of the attractive interaction. Thus, as demonstrated for three dimensional 59,60 and two dimensional colloidal systems dispersed in simple liquids, [61][62][63][64][65] using computer simulations and experiments, we expect that NIs aggregate is of fractal dimension.
In Fig. 3, we report RDF versus the dimensionless distance, r* ¼ r/s, for several values of the NIs density, using MD simulation. Firstly, we remark that at a small distance, compared to the NIs diameter that is for 0 < r < s, RDF practically vanishes due to the imposed strong steric repulsion between the adjacent NIs, which is applied to prevent the nanoparticles from overlapping. Then, for 0 < r < s, the NIs present a correlation hole. For distances beyond s, the NIs show a local order. The latter is due to a strong correlation between NIs. Of course, the size of this region crucially depends on the value of the physical parameters appearing in the interaction potential, in particular, the density of NIs. At a very high distance, however, RDF saturates naturally to 1, whatever is the value of density. This limit means that the NIs are completely disordered in space, i.e., the order is absent at a long-range. Between the two distance-regimes, RDF oscillates and presents as a succession of peaks of decreasing heights that indicating the existence of many shells of neighbours. The locations of maxima of these oscillations correspond to the preferred relative distances between the neighbouring NIs. Secondly, as shown in Fig. 3, the position of the RDF rst peak, r max , is shied toward the le as the density of NIs increases. This shi is due to a net dominance of the attractive part of the mean-eld interaction potential, which is related to the RDF by g(r) ¼ exp(ÀU/k B T), as shown in Fig. 4. Then, the coordination number, N c , increases progressively, as observed in Fig. 5.
From Fig. 3, we extract the reduced maximum position of the principal peak, r max /R, for a wide range of NIs number density. A good t with the calculated numerical data from MD simulation shows that r max decreases with the NIs number density according to a straight line of the equation, r max =s ¼ À937:7n þ 4:8; n # n c r max =s ¼ À29:45n þ 1:58; n $ n c Also, Fig. 3 shows that the principal peak-height increased by an increase of the NIs density. We attribute such an increase to a progressive dominance of the attractive part of the mean-eld interaction potential of Fig. 4. Thus, the latter gives rise to an aggregation of the NIs in the host membrane. The increase of the maximum, g max , as a function of r max /s is depicted in Fig. 6. The best t of numerical data gives the following form for g max .
lnðg max Þ ¼ À0:30 lnðr max Þ þ 0:46; n* # n c lnðg max Þ ¼ À10:47 lnðr max Þ þ 4:5; n* $ n c The MD calculation of d F is in good agreement with the Scaling Theory, it suggests that for number densities below the gelation one n < n c , d F does not depend on n. In this sense, it is a universal quantity. However, as we will see later, the Scaling theory fails to predict the d F for number densities above the gelation density n > n c , giving negative value of d F . We add that the amplitude, H, can be obtained by extrapolation on the vertical axis. A comparison with the scaling-relation means that the NIs aggregate is fractal and we have, for n # n c , d F ¼ 2 À 0.30 ¼ 1.70. This value is slightly great than the fractal dimension calculated for simple colloidal systems, d F ¼ 5/3, 63 that found in good agreement with other experimental results. 62,65 These results inspire a more general understanding of the liquid-gel transition in colloidal systems. To discuss the possibility of a reverse gelation process, we recall that the production of amorphous solids from an ordered solid has been identied. 66,67 In practice, this transition can be made under the effect of heating or stretching the materials, 68,69 which will increase the volume and therefore decrease the density. The ow process, which thins an initial solid material, can be seen as a process opposite to gelation, which solidies an initially liquid sample. It seems that the observed appearance of stiffness from a uid state and the onset of the ow of a solid-state are two almost mirroring manifestations of the phase transition. 21

Thermodynamic properties
The aim of this section, recalling that we conducted the simulations in xed temperature, is to study the effect of NIs number density on the microscopic potential energy, E pot , and the pressure, Press.
We calculated the potential energy as: where U eff (ij) is the effective interaction potential dened above. We calculated the pressure as: where the rst term represents the kinetic pressure due to the NIs choc with the walls of the simulation square, like an ideal gaze, and the second term is the Virial Pressure, which is due to the inter-particles choc. Fig. 7a shows the variation of the potential energy (E pot ) as a function of the NIs number density (n). We note that the E pot is close to zero for very low densities but shows a logistic growth when (n) increases. For n # n g , E pot shows a rapid linear increase. For n > n g , E pot shows a slow increase. But when the numerical density reaches the gelation density (n ¼ n g ), E pot shows a sharp jump. The increase in E pot is due to an increase in the correlation between particles, as quantitatively shown in Fig. 3-5. Fig. 7b represents the variation of the pressure (Press) as a function of n. The pressure gradually increases with the numerical density, with a relatively faster increase for n > n g , and as well as E pot , shows a noticeable jump for n ¼ n c . The numeric values of E pot and Press are given in Table 1.

Dynamic properties
3.3.1 Theoretical prediction of dynamic properties for low density nano-inclusions. Nano-inclusions experience a ballistic diffusion, for early times, because of the absence of the correlations. At these small-time-scales, the random walker is not yet gets stuck in a cage. As we shall describe below, such a regime is put in evidence by our MD simulations, at the short times.
We describe the Brownian motion of a NIs by the following phenomenological Langevin equation: where m stands for the mass of the tracer, v(t) for its velocity, z for the friction coefficient, and F s (t) for a random force. The friction coefficient z is related to the viscosity of the NIs, h and the NIs radius, R, by the classical Stokes relation z ¼ 6phR. Of course, the viscosity of the solution depends on the number density of the dispersed nanoparticles n. For small values of the number density, this effective viscosity obeys the well-known law, 70 here n accounts for the reduced NIs number density, h M is the viscosity of the membrane, and h N is that of NIs. As it should be, we have h(n) > h M , since h N > h M . We return to the stochastic force, and we see that we considered as white noise with: where d(t) is the Dirac distribution. The brackets h.i mean an average over time.
In the following, we are interested in the mean-squaredisplacement (MSD) dened as h(r À r 0 ) 2 i h hD 2 ri, such as: with g ¼ z/m knowing as the relaxation rate. For times much longer than the inverse relaxation rate, that is for t [ g À1 h s, MSD grows linearly with time, and depend on the diffusion coefficient D, dened as the compromise between the thermal uctuations and the dissipation such as D ¼ k B T z in two dimensions, by:  For t ( s the Taylor development of e Àgt gives: e Àgt z 1 À gt þ 1 2 ðgtÞ 2 , then we nd: 3.3.2 Theoretical prediction of dynamic properties for high density nano-inclusions. When the number density of NIs is riched n c their dynamics are more complex than those in low number density. In this case, NIs particles move ballistically at short times. Thus the mean-square displacement hD 2 ri $ t 2 , which is followed by a crossover to Fickian diffusion, characterized by hD 2 ri $ t for long times. But, in dense colloidal systems a caging effect, where the atoms are trapped by their neighbors, cause a subdiffusion regime that intermediate the ballistic and the diffusive regimes. This regime is characterized by hD 2 ri $ t a , with 0 # a # 1. a goes to a limit value when the number density reaches the gelation one that we note n g . The analytical determination of the value of a for the subdiffusive regime remains uncertain. However, to predict an approached value of this key parameter, there are developed theoretical approaches based on the inclusion of complex phenomenological memory functions in the generalized Langevin equation. In this context, Zwanzig has developed a Dynamics Theory based on the generalized Langevin equation (GLE) instead of the standard Langevin equation (SLE). 71,72 This approach can be adopted to discuss the subdiffusion observed in NIs, which we considered as two-dimensional colloidal particles swimming in a solvent (the host membrane). From a mathematical point of view, it is an extension of the SLE developed to study the dynamics behaviour of simple liquids, where the friction is assumed to be determined using the instantaneous velocity of the particles under study. We note that the difficulty lies in understanding the phenomena of subdiffusion is to explain the cage effect mathematically. Besides, physically speaking, the subdiffusion depends on various parameters, namely, temperature, number density, and nally, and least lost the nature of the host medium. In this context, we express the GLE as, where v the velocity of a moving nano-inclusion (a tracer).
g ¼ z m is the relaxation rate, where z is the friction coefficient, and m is the mass of the tracer. k is the memory-function that expresses the friction retardation, and F s (t) is a random force felt by the moving nano-inclusion due to its collisions with the molecules of the host membrane. Thus, the random force verify the equation, For these considerations, VACF solves the following differential equation: For the analytical resolution of the GLE equation, we refer to an elegant memory function, proposed by Flenner et al., 73 originally developed to study the dynamics aspect of lipid  atoms, which also presents a subdiffusion phenomenon. This approach base on the Zwanzig-Mori projection method to model the hD 2 ri behaviour. This is done by involving the equation of motion for the density autocorrelation function F s (q, t) ¼ hn(Àq, 0)n(q, t)i of a tracer, i.e. a selected particle, at the wave vector q, where v ¼ ffiffiffiffiffiffiffiffi k B T p =m, k B stand for the Boltzmann's constant, T denote the temperature and m refer to the nano-inclusion mass.
We obtain the two-dimensional equation of motion of the MSD from: where k 0 ðtÞ ¼ lim Mðq; tÞ.
The Flenner et al. approach consists of introducing an adequate memory function that explicit three crossover times scales and thus reproduce the ballistic, the normal, and the subdiffusive regimes, 73 where d(t) is the Dirac delta function, B is a dimensionless parameter. s 1 is the characteristic time for the crossover from the subdiffusion to the normal diffusion. s 2 is the onset time of the subdiffusion regime. And, s 3 is a characteristic time for the crossover from the ballistic to the subdiffusive regime. We note that, for B ¼ 0, we recover the standard Brownian diffusion laws, and the hD 2 ri behaves like: If (t/s 2 ) a ( 1, then The power-law term of the memory function represents the subdiffusion phenomenon. Flenner et al. performed a numerical computation, taking into account this term, 73 and revealed that there is a subdiffusive regime with a < 1 intermediating the ballistic and normal regimes. This region is due to a cage effect applied to a considered particle by its surrounding neighbors. Thus, the mean-square-displacement, hD 2 ri, is expressed as: 3.3.3 Analysis of MD simulation results. In the last years, MD simulation becomes a privileged method for the statistical description of dynamic properties. In the framework of this method, the dynamics behavior is characterized quantitatively by the mean-square-displacement (MSD), the velocity autocorrelation-function (VACF), and the diffusion coefficient (D a ).
In this section, we shall discuss the inuence of the number density variation on the dynamic properties of the BuSn12 membrane NIs, using NST-MD simulations. We recall that the considered NIs are characterized by a radius R ¼ 4.5Å and molecular weight of 2866.7, the temperature xes to room temperature, but the number density of NIs varies on a wide range.
In Fig. 8, we present the log-log curve of MSD as a function of reduced time, for different nano-inclusion number density values. For low densities (n < n c ¼ 3.67 Â 10 À3ÅÀ2 ), we note that the MSD have two diffusion regimes characterized by a transition time (t # s), i.e. a ballistic diffusion with an exponent (a ¼ 2) for t # s and a normal diffusion (a ¼ 1) for t > s, we notice that s decreases when n increases, that is to say that the diffusion of the nano-inclusion quickly reaches the normal diffusion. Moreover, for high densities n > n c , we note that the MSDs curves have three diffusion regimes. At more important times, there is a crossing towards the Fickian diffusion with hD 2 ri ¼ 4D a¼1 t, where D a¼1 is the long time self-diffusion coefficient. But, in this crossing region, from the ballistic to the normal diffusion, a subdiffusive region appears where hD 2 ri ¼ 4D a t a , with a < 1. The value of D a<1 is smaller for high density and has a limit value of about 0.5 for the rst gelation number density (n g ¼ 8.40 Â 10 À3ÅÀ2 ). Fig. 9 describes the variation of VACF as a function of time for several number densities, below and above n c . We characterize this function by more pronounced damping oscillations at high density. This behavior is in good agreement with the results of the theory described above. In high density, we attribute the VACF oscillations, which take negative values, to the connement of each nano-inclusion in a cage formed by their close neighbors. These oscillations signify that the diffusion has a memory. On the other hand, VACF is positive and has few damping oscillations for low densities as a comparison to the high-density systems. Also, we note that the asymptotic behavior of VACF, which reects the diffusion regime of a studied system, tends towards zero negatively for high-density values, according to the theory described above.
In Table 2, we give the diffusion exponent, a, and the corresponding crossovers time for n < n c . We note that a takes two values, 2 and 1, corresponding to the ballistic and normal regimes, respectively. The crossover time from the ballistic regime to the normal one decrease as the number density increase. This decrease is due to an increase in the correlation between NIs.
In Table 3, we give the diffusion exponent, a, and the corresponding crossovers time for n $ n c . We note that a takes three values, namely, the two values observed for n < n c (2 and 1), in addition to the value a < 1 corresponding to a subdiffusive regime, which intermediates the ballistic and normal regimes, this is due to the cage effect discussed above. We note that, for the subdiffusive regime, a is not unique for all densities. In fact, a decrease progressively and reaches a limit value a ¼ 0.5 for n g .
In Table 4, we give the normal diffusion coefficient of BuSn12 NIs, D a¼1 , in reduced units and SI units, as a function of n. In Table 5, we give the subdiffusion coefficient D a<1 , corresponding to the subdiffusive regime. We note that the subdiffusion appears for a small time-interval except for the gelation number density n g the NIs dynamics is blocked for a long time. Fig. 10a and b depict the diffusion coefficients D a¼1 and D a<1 , for the BuSn12 NIs, giving in Tables 4 and 5, respectively.
We report that the variation of the D a¼1 and D a<1 as a function of the number density n follow a rst-order exponential decay (Arrhenius-like decay).
For D a¼1 we have, with the best tting parameters D 0 ¼ 4.16 Â 10 À9 m 2 s À1 AE 1.67 Â 10 À9 m 2 s À1 (39) and for D a<1 we have, with the best tting parameters D 0 ¼ 7.96 Â 10 À11 m 2 s À1 AE 1.77 Â 10 À10 m 2 s À1 (43) We note certain limitations of this study. First, although the calculated properties follow well-established physical laws, the tted parameters are not universal and depend on other relevant factors, such as the shape and size of the NIs. Therefore, projecting the results of the present study on non-spherical NIs may be inappropriate. Second, we also planned to explore the physical properties at n > n g , but the MD simulation code is proving to be unstable. Possibly the MD divergence is due to the attractive part of the interaction potential.
From the literature, for colloidal suspensions in simple solvent (water, for example), it is known that the subdiffusion phenomenon occurs at high enough number density, which induces the cage effect, 75,76 where the diffusion of a selected particle getting trapped in a cage formed by its close neighbors. In contrast, in the present case, we mark that this phenomenon appears for low number densities. It is natural to link the membrane NIs slow diffusion to the nature of the host medium. The latter is a biological membrane composed of highly ordered lipids that tend to limit their diffusivity.
To compare the diffusion of proteins and the diffusion of inorganic NIs within biological membranes, we refer to our recent work. 36 In this previous work, we have studied the subdiffusion of cylindrical shaped proteins inclusion in phospholipids mixture bilayer, using coarse-grained MD simulations. Accordingly, for the room temperature, we found a mean subdiffusion exponent of proteins to be a ¼ 0.65, and we found the diffusion coefficient of the order of 10 À16 m 2 s À1 this range is also found in some experimental works. [77][78][79] In the present work, the diffusion of inorganic nanoinclusion is in the range of 10 À14 to 10 À9 m 2 s À1 . These values overlap with those of other similar works, which is in the order of 10 À14 m 2 s À1 , 80,81 performed for a number density close to n g .
Thus, we deduce that the diffusion of inorganic NIs is faster in comparison to the organic inclusions.
We recall that synthetic membranes can be composed of single lipids type, the liposome membranes as an example. 82,83 But the biological membranes are crowded structures encompassing several lipids types, included and external proteins, and present in different phases, the pulmonary surfactant as an example. 36 Thus, the expensive details related to the membrane nature limit the studied number of NIs in an MD simulation. For this purpose, a direct interaction potential between NIs, which respect the membrane uctuation, permits the study of a large NIs number. Thus, an effective interaction potential cannot provide nanoscopic properties. Although, it gives an overview of the NIs mesoscopic properties. In this context, the NIs are considered a two-dimensional colloidal dispersion system. For simplicity, the nano-inclusions represent the colloidal dispersion, and the membrane is the two-dimensional solvent.
In this work, we studied the BuSn12 nanoparticles as membrane nano-inclusions, and we assumed that they are hard-spheres and charge neutral. We perceive that their interaction-potential form with the appropriate parameters is  determined experimentally in previous works by D. Constantin. 37 Generally, nano-inclusions can be charged or neutral and present in different shapes and sizes, depending on their applications. [84][85][86][87] From a mesoscopic viewpoint, MD simulations require appropriate interaction potential forms and parameterizations. These proposed models are oen semiempirical. For instance, the electrostatic interaction due to the surface charge can be modeled using a Yukawa screened coulombic interaction type, mimicking the standard DLVO model for charged colloidal particles in a simple solvent. 76,88 Besides, the membrane uctuation can also depend on the nano-inclusions shape. Thus, the interaction potential form and parameters are specic to each nano-inclusion type and cannot be precisely a Gaussian type.
Consequently, the three scientic approaches encompassing theory, experiment, and simulation are needed to explore the biophysical-chemistry prospects of various nano-inclusions types at different time-length scales.

Conclusion
In this work, we performed MD simulations in the NST statistical ensemble conditions to study membrane NIs, in which a specic type of nano-inclusion, BuSn12, is discussed in detail. We conducted the study for a wide range of NIs densities. As hypothesized, the calculated physical properties have shown a remarkable discontinuity, corresponding to the critical and gelation number densities. Thus we have provided a mathematical understanding of the spherical NIs aggregation.
While previous works study the interactions between BuSn12 NIs across the membrane, this study offers further physical insights into their phase transition phenomena. Accordingly, we discussed the structural properties relative to the number density by developing a Scaling Theory approach. This approach base on the idea that the NIs aggregate is of fractal structure. Such a scaling method consists with the MD simulation for n < n c but diverge for n > n c . We determined the thermodynamic properties through the potential energy per nano-inclusion and the pressure, using MD simulation. Thus, we observe that these major physical quantities present a discontinuity for n ¼ n c and a clear jump for n ¼ n g . We studied the dynamic properties and their dependence on the NIs number density through the temporal evolution of the mean-square-displacement, the velocity auto-correlation-function, and the diffusion coefficients. Precisely, we explore the cage effect and the subdiffusion phenomenon in membrane NIs, using two approaches, which are a Dynamics Theory, based on the generalized Langevin equation, with the Flenner memory function, and MD simulation. We observe that the two dynamics approaches are in good agreement. Thus, the thermodynamics and dynamic properties analysis allow determining n c and n g .
Gelation is a common phenomenon in colloidal science, which has several nanotechnological and medical implications. Thus, the simulation protocol we performed to study spherical membrane NIs, using a solvent-free interaction potential, can be adapted to different shapes and types of NIs, such as nano-rods and nano-discs. The approaches we used to determine the phase transition parameters can be conducted for other colloidal systems, such as microemulsions and metallic nanoparticles. In particular, to reveal drug delivery and cancer treatment.

Conflicts of interest
There are no conicts to declare.