Monte Carlo Simulation of Heterotypic Cell Aggregation in Nonlinear Shear Flow

In this paper, we develop a population balance model for cell ag-gregation and adhesion process in a nonuniform shear flow. Some Monte Carlo simulation results based on the model are presented for the heterotypic cell-cell collision and adhesion to a substrate under dynamic shear forces. In particular, we focus on leukocyte (PMN)-melanoma cell emboli formation and subsequent tethering to the vascular endothelium (EC) as a result of cell-cell aggregation. The simulation results are compared with the results of experimental measurement. Discussions are made on how we could further improve the accuracy of the population balance type modelling. 1. Introduction. Neutrophil and tumor cell adhesion in a nonuniform shear flow is a problem of great interest in the study of tumor extravasation during metastasis. There have been many computational and experimental studies of cell aggregation in the flow environment. Detailed investigations of the aggregation and adhesion process involving individual cells have been done for both uniform and nonuniform Statistical population studies have also been done for the constant shear flow case. For example, Laurenzi and Diamond (see [3]) used the Monte Carlo method to simulate platelet and neutrophil heterotypic aggregation in a cone-plate linear-shear assay. Parallel plate flow chamber experiments have shown that under different flow conditions, melanoma cells and leukocytes (PMN) adhere to each other with different efficiencies. However, compared with the case of constant shear rate, the study of aggregation of cells in nonlinear shear flows poses an even greater challenge. We have undertaken a study to explore how the interactions between PMNs and melanoma cells are affected by the fluid dynamics in nonuniform shear flows. A special focus is on the PMN melanoma cell emboli formation in a nonuniform shear flow and subsequent tethering to the vascular endothelium (EC) as a result of cell cell aggregation. As the first step, we developed some computational tools which can be used to conduct statistical studies of heterotypic cell cell collision and adhesion to


Introduction. Neutrophil and tumor cell adhesion in a nonuniform shear flow is a problem of great interest in the study of tumor extravasation during metastasis.
There have been many computational and experimental studies of cell aggregation in the flow environment.Detailed investigations of the aggregation and adhesion process involving individual cells have been done for both uniform and nonuniform shear flows in [1,5,6,7,18,19,24,25,29].Statistical population studies have also been done for the constant shear flow case.For example, Laurenzi and Diamond (see [3]) used the Monte Carlo method to simulate platelet and neutrophil heterotypic aggregation in a cone-plate linear-shear assay.Parallel plate flow chamber experiments have shown that under different flow conditions, melanoma cells and leukocytes (PMN) adhere to each other with different efficiencies.However, compared with the case of constant shear rate, the study of aggregation of cells in nonlinear shear flows poses an even greater challenge.
We have undertaken a study to explore how the interactions between PMNs and melanoma cells are affected by the fluid dynamics in nonuniform shear flows.A special focus is on the PMN melanoma cell emboli formation in a nonuniform shear flow and subsequent tethering to the vascular endothelium (EC) as a result of cell cell aggregation.As the first step, we developed some computational tools which can be used to conduct statistical studies of heterotypic cell cell collision and adhesion to a substrate under dynamic shear forces, along with in vitro experiments.The main objective of this paper is to report our progress in the application of population balance models and to discuss the need for further analysis and modifications.
In our in vitro experiment [20], the fluid flow in the parallel plate flow chamber assays can be described as Poiseuille's flow through a rectangular geometry, which characteristically has a parabolic velocity profile, as shown in Figure 2, instead of a linear shear rate.In the case of melanoma cells and PMNs collision/aggregation, a vessel wall or experimental substrate is always present and therefore changes the hydrodynamics of the system, and thus the collision probability (see Figure 1, courtesy of the Cellular Biomechanics Laboratory at PSU).We are interested in developing a mathematical model to simulate cell aggregation in the near wall region.At the microscopic level where the collision and adhesion of two cells takes place, the aggregation depends mainly on two parameters, namely, the collision rate and the adhesion efficiency.The collision rate characterizes the probability that two given particles will collide, which may depend on the shear rate, the diffusion, the gravity force and the particle size.Two cells adhere due to receptor ligand bonds.The adhesion efficiency is the probability that two cells will adhere to one another when they collide, and it reveals how this probability may be affected by numerous factors, such as the shear rate, the particle size and the particle surface properties, etc.Our aim is to use numerical simulations based on population balance models as a tool to investigate the roles played by the nonuniform shear field, the collision rate and the adhesion efficiency in the aggregation process, compared with the in vitro experiments.The main mathematical model used in this study is the population balance equations (PBE), or the Smoluchowski equations.The details are presented in Section 2. Numerical simulations based on the PBE provide statistical information on the cell aggregation and adhesion process.The detailed description of the simulation methodologies are presented in Section 3, which is the main focus of this paper.One advantage of numerical simulations based on the population balance type of models is that the much of our numerical results can be readily compared with the in vitro experimental observations.This is useful for parameter calibration, model validation and model prediction.We also derive some interesting scaling invariance properties of the population balance equation, so that they can be utilized to better fit the experimental conditions.
2. Population balance model.Shear induced aggregation is a phenomenon that, in the case of platelets or PMNs, has been analytically and numerically modelled [12,15,16,17,22,26].Of particular interest to us is the use of statistical type models for predicting and simulating cell aggregations, namely, the population balance models.In such models, the aggregation of two individual cells mainly depends on two parameters: the collision rate and the adhesion efficiency.The collision rate characterizes the probability of two given particles colliding, which is a function of shear rate, particle size and convection.The adhesion efficiency is a measure of the probability of two cells adhering as a response to the collision, which is a function of shear rate, receptor/ligand density, and receptor/ligand avidity.The population balance equations, as introduced by Smoluchowski [2], are also referred to as the coagulation, or Smoluchowski, equations and have been applied to a wide range of applications such as aerosol growth, polymerization problems, and the kinetics of platelet aggregate formation and disaggregation.
In mathematical terms, the population balance model describes a rate of change of the density of a particle of a particular size as a function of time: where C(x, t) is the concentration of the particles of size x at time t, and β(x, y) is the coagulation kernel which describes the coagulation probability between two particles with size x and y.The first term of the right hand side describes the generation of the particle of size x by the aggregations of smaller particles.For example, if a particle of size x − y adhere to the particle of size y, it becomes a particle of size x.The second term describes the loss by aggregation with other particles.
The coagulation kernel, originally derived for modeling collisions in laminar shear, contains a constant shear rate γ [2].The basic idea is to consider a moving particle of radius x sticking to the one of radius y in the linear shear flow.Then we compute the number of point masses which hit the sphere of radius x + y per unit time.The computation shows that the collision kernel has the following form: Such a kernel was designed to model systems that did not contain cells; thus modifications of the original coalescence kernel have been proposed in the literature to accommodate kinetics properties of receptor ligand type binding in cellular systems.Generally, a term referred to as the adhesion efficiency ε, has been introduced, resulting in a new kernel of the following form (see [12,15,16,17,22,26]): where ε is the adhesion efficiency and γ is the shear rate of the linear shear flow.The adhesion efficiency term was first estimated by Belvel and Hellums [12] and later extensively studied by Huang and Hellums [15,16,17].In their model, key intrinsic biological parameters were estimated by matching the theoretical results with an experimental volume distribution curve of shear induced platelet aggregates.The experimental data was obtained by shearing platelet suspensions in a cone plate viscometer and analyzing the size distribution of the aggregates by Coulter counter.In later works [22,26,27], similar methods were also used to estimate adhesion efficiency by fitting experimental data to equations based on the Smoluchowski coagulation theory.In such analysis, the estimated adhesion efficiency term was theoretically de-convoluted into a receptor component and a hydrodynamic component.Yet, the existing analyses have focused on the cell collision/aggregation in a linear shear flow far away from a boundary by using cone plate viscometer assays.However, in our experiments, due to the nonlinear shear flow condition, the aggregation kernel that they used for linear shear flow is not appropriate anymore.The adhesion efficiency will be different under different flow conditions.Thus, such approaches are not appropriate to model the experimental setups that are considered here, i.e., in characterizing cell cell interactions in the near wall region, which are more important in understanding the overall behavior of adhesion mediated melanoma extravasation.Our task is to look for a more suitable aggregation kernel, measure the adhesion efficency in the near wall region, etc.
Other researchers have derived alternative coalescence kernels to accommodate different hydrodynamic conditions, for example, in turbulence and gravitational settling.One particularly applicable form of the kernel was developed to model particle collisions in a nonlinear shear field [11].It was derived by considering the collision frequency of droplets in a turbulent flow field [28] to have the following form: ) where the shear rate and the adhesion efficiency are not included in this equation.However, these and some other constant parameters have been lumped into a single term, c.
Some existing studies have suggested PMN mediated melanoma cell extravasation is influenced by cell populations in a shear flow.In order to understand how heterotypic cell ratios affect melanoma PMN collisions and subsequent aggregation in the near wall region (hence affecting PMN facilitated melanoma extravasation) under flow conditions, we have begun to work with a special variation of the population balance (PB) model, as introduced by Smoluchowski.By comparing the simulations to experimental results, the model can be validated.Through both ad hoc and systematic modifications to existing models, the PB equation can accommodate changes in collision rate and adhesion efficiency due to changes in hydrodynamic parameters.The variation of the PB model studied here allows the simulation of the heterotypic cell cell aggregation in a nonuniform shear field, especially the aggregation in the near wall region.This provides a potentially useful method to model the PMN melanoma collision and aggregation under a wider range of more physical conditions than those that can be reasonably tested experimentally.A more physiologically relevant case will be when the number of PMNs is significantly greater than melanoma cells (normal physiological levels are 2 × 10 9 − 7.5 × 10 9 PMNs per liter of human blood).We will revisit this issue later in the discussion section.

3.
Simulations of the PB model, methods and parameters.
3.1.Heterotypic aggregation near the wall.As discussed previously, most existing experimental and mathematical applications of PB models for analyzing cell cell collision/aggregation have been based on spatially homogeneous coagulation in a linear shear flow.Other researchers have derived alternate forms of the coagulation kernel to accommodate particle collisions in a nonlinear shear field [11], as shown in equation ( 4).An earlier study published by Laurenzi and Diamond [3] provided an example of modeling heterotypic cell cell collisions.In this model, shear induced PMN platelet aggregation was predicted by the PB equation, which was solved by a Monte Carlo (MC) driven algorithm.Similar to the Laurenzi and Diamond study, the input parameters for the melonoma PMN collision model were chosen to be cell density, cell size, wall shear rate and approximate adhesion efficiency.The output was the percentage of aggregate formation at a given time.
There are two main problems that need to be solved in developing our PB model in order to match our experiments.First, due to the nonlinear shear flow, the concentration in the near wall region will be different from the inlet cell concentration in the chamber.However, we can only use the inlet cell concentration as our model initial data since it is very difficult to measure the concentration in the near wall region in the experiments.Fortunately, we can solve this problem through scaling invariance.
The second problem is to modify the existing coagulation kernels such that it can describe the collision in nonlinear shear field, especially in the near wall region.These two problems will be discussed in the following subsections.
3.2.Scaling invariance.As we indicated, we will use the cell concentration in the inlet of the chamber as our initial concentration in our numerical simulation algorithm, since it is very difficult to measure the concentration in the near wall region in the experiments.It is evident that the local concentration in the near wall region will not be the same as the inlet concentration due to the flow action; thus, we need to scale it such that these two are consistent.Here we make a simple assumption, that is, that the ratio of the concentration in the near wall region to the concentration in the inlet of the chamber is a constant.Then, taking this into account, we may normalize our initial data mathematically by using the scaling invariance of the population equations.
The scaling invariance refers to the following mathematical properties.Let C(x, t) be the solution corresponding to the initial concentration C 0 (x) and λ be any positive parameter; then we define a new function C(x, t) = λC(x, λt). ( One can easily see that C(x, 0) = λC 0 (x).Moreover, direct calculation gives This shows that C(x, t) = λC(x, λt) is a solution of the same population balance equation with a scaled initial condition.This is an interesting scaling invariance property of the population balance model.An alternative way to represent the scaling invariance without rescaling time is to explore the normalization of the coagulation kernel.We can also verify that if the coagulation kernel β(x, y) is multiplied by a positive constant λ, denoted by β(x, y), then the function gives a solution of the population balance model corresponding to the kernel β(x, y) with a scaled initial condition λ −1 C 0 (x).
Of course, based on the above discussion, one can also obtain the scaling invariance by rescaling both the kernel β and the time while keeping the same initial condition.
An important conclusion based on the above observation is that we can modify the kernel function by a constant factor in order to renormalize the total initial concentration.It should be pointed out that the non uniform flow conditions not only alter the total concentration of the cells near the wall, but also have different effects on different types of cells.The simple renormalization of the coagulation kernel cannot account for this further complication, which leads to an issue discussed later.Here we will use the kernel which has the form similar to the one given by equation ( 4).But we will combine the shear rate and adhesion efficiency into this kernel.To estimate the shear rate that melanoma cells and PMNs likely experience in a near wall region, an average shear rate was calculated.This average shear rate was determined to be the average shear rate from one cell radius (the closest a cell can be to the wall without penetrating the wall) to four cell radii from the wall (two cell diameters) (Fig. 2).Note that nonlinear shear flow has a quadratic velocity distribution in the flow chamber [29]: where y is the y-axis coordinate, Q is volumetric flow rate, b is half the chamber height and w is flow chamber width.Thus, the shear rate along the y-axis has a linear distribution: Then we get the average shear rate in the near wall region: where a is cell radius.Then, utilizing this local shear in the coagulation kernel, we have: where we have introduced a renormalization function φ, which is a function of the local melonoma and PMN concentrations, C M and C P M N , local cell Reynolds R e and Capillary numbers (deformation)C a and wall proximity, δ.The renormalization constant will be discussed in sections 3.5.

3.4.
Determining the adhesion efficiency.An important parameter used in equation ( 8) is the adhesion efficiency, ε.Heterotypic collision models require three separate values.In the case of melanoma cells and PMN collisions, the relevant ε are ε T +P , ε P +P and ε T +T .The most general definition for adhesion efficiency is the number of formed aggregates divided by the number of collisions that occur.In the model developed here, the values for PMN PMN (P + P ) homotypic adhesion efficiency, ε P +P , were taken from the same published source as Laurenzi and Diamond [3] (Table 3.4), and the value for ε T +T was approximated to be zero based on our own experimental observation [21].
To determine melanoma PMN (T + P ) adhesion efficiency, ε T +P , aggregation data from a parallel plate flow assay was used (see Table 3.4).Adhesion efficiency was calculated by dividing the number of melanoma PMN aggregates formed by the number of melanoma PMN collisions that were counted within the same field of view.The adhesion efficiencies were determined for three cases: untreated, anti-CD11a and anti-CD11b.CD11a and CD11b are two different α chains, involved in β 2 integrin expressed on PMNs.CD11b/CD18 is called LFA-1 and CD11b/CD18 is called MAC-1.Table (3.4)shows the ε P +P and ε T +P derived from experiments.1), assuming a spatially homogenous 1 : 1 ratio of C T C : C P M N , using the new kernel (8), and running the simulation without normalizing the initial population.The predicted percentage of melanoma PMN aggregation near the substrate was then compared to the experimental percentage of aggregation determined from the parallel plate assays.The normalization factor φ, necessary to correct the initial concentration in the near wall region to match the experimental results, was then calculated.Its value was derived using the untreated melanoma and PMN cases as a control data set.The parallel plate experimental percent of aggregation was determined by the number of melanoma cells that adhered to the substrate due to a melanoma PMN collision divided by the flux of melanoma cells near the surface [21].This normalization by the flux of melanoma cells near the surface was necessary to account for the differences in cell flux with varied shear rate [23].
3.6.The numerical algorithm: a Monte Carlo method.Some numerical methods have been developed to solve the population balance equation, we refer to [14] for additional references.Here we use the Monte Carlo method introduced by Gillespie [8,9,10].
Since we consider the heterotypic cell aggregation, the kinetics of multi component aggregation process will be described by the discrete population balance equation as follows: β(i, j; i , j )C(i, j; t)C(i , j ; t). ( Note that all the properties of the continuous PBEs we discussed previously are all true for the discrete PBEs. In our case, C(i, j; t) denotes the concentration of particles with i melanoma cells and j PMNs.β(i, j; i , j ) is the coagulation kernel corresponding to the adhesion event A i B j + A i B j → A i+i B j+j , where A and B represent the melanoma cell and PMN, respectively.
We will denote the composition of particle A i B j by [i, j] later on.Gillespie's algorithm (see [8]) reveals that the coagulation frequency γ(i, j; i , j ) in the population balance equation has the property that γ(i, j; i , j )δt = probability that a given pair of particles with compositions [i, j] and [i , j ] will aggregate in the next time interval δt.
The relationship between the coagulation frequency and the coagulation kernel is established in [8]: where V is the total volume of the aggregation system.Consider the aggregation between two particles with compositions [i, j] and [i , j ] in an aggregation system.If the aggregation event is the unlike interaction, then the aggregation probability in the next time interval δt is a([i, j]; [i , j ])δt, where If the aggregation is the like interaction between the two particles of one species where X [i,j] and X [i ,j ] are the populations of species [i, j] and [i , j ], respectively.For the Monte Carlo stochastic approach, the key point is the transition from one state of the aggregation system to the next one.The aggregation probability density function, P ([i, j]; [i, j]; t), is defined as the probability that two particles of species [i, j] and [i , j ] will aggregate in the next time interval δt after an interval of quiescence (0, t).
Let P 0 = P 0 (t) be the density function for the imminent quiescence time t.By definition, P ([i, j]; [i , j ]; t)δt is the probability density function that there is no aggregation in the interval [0, t], with the aggregation event only taking place in ) .Also note that for a small time interval [t, t + δt], the probability P 0 (t + δt) that the aggregation system remains in the state will be is the probability that there is no aggregation event in [t, t + δt].Therefore we get an equation for P 0 (t) as follows. where Solving this equation, we have P 0 (t) = e −αt .Then we get the probability density function, as derived already in [4,8,9,10], We may rewrite the aggregation probability density function as: where α is the total aggregation probability frequency defined by equation ( 14).Thus we can get the probability of a particular type of aggregation in the time interval t.
) By the Monte Carlo method, we generate two random numbers, r 1 and r 2 uniformly in the interval [0, 1], to calculate t, the aggregation period, and the aggregation event between species µ and ν.Here t is generated by equation (18), The aggregation event between species [i, j] and [i , j ] is generated in the following way.Assume the aggregation system has m species at a certain time.Let us label the species by the number from 1 to m.Thus for any two species [i, j] and [i , j ]), they correspond to a pair of numbers, (p , q ), respectively.Here the aggregation event between two species, (p, q), is determined by inequality (19), Therefore the aggregation event will be determined.These calculations lead to the following simulation steps.

Algorithm:
1. Given N initial species and their populations, P i , calculate their aggregation probabilities by the aggregation kernel β(i, j).Compute the total aggregation probability, α. 2. Generate two random numbers, r 1 and r 2 , determine the aggregation time step, t, by equation (18), and calculate the aggregation event by (19).3.After the aggregation takes place, update the number of species and the populations.4. Go back to step 1 until the aggregation stops.The consistency of the Monte Carlo algorithm with the PBE models can be easily established.In the population balance equation, the term β(i, j; i , j )δt/V represents the coagulation probability between two given particles with components of [i, j] and [i , j ], respectively, in the next time interval (t, t+δt).This is consistent with equations ( 11) and (12), which represent the coagulation probability between two species.We will show that the population balance equation can be derived from the stochastic aggregation model contained in this Monte Carlo method.
For an arbitrary given particle, we assume its composition at time t is [i, j].Thus according to equation (10), the probability that this particle will adhere to any particle with component [i , j ] in the next time interval (t, t + δt) should be β(i, j; i , j )X(i , j ; t)δt/V , where X(i , j ; t) is the population of the particles with component [i , j ] at time t.The probability that it will not adhere to any particle is thus 1 Assume the density distribution function of the particle with component [i, j] at time t is f (i, j; t), which is defined as where v 1 and v 2 are the volume of tumor and PMN.Then it satisfies β(i, j; i , j )X(i , j ; t)δt   f (i, j; t) β(i − i , j − j ; i , j )δtX(i , j ; t)f (i − i , j − j ; t) .

Thus we have
β(i, j; i , j )X(i , j ; t)f (i, j; t) Take the limit of the left hand side when δt goes to 0, and we have β(i, j; i , j )X(i , j ; t)f (i, j; t) Note that for any two constants v 1 and v 2 , we change the index i − i = ĩ and j − j = j; then we have

So we get
Note that the concentration C(i, j; t) is defined as Then, dividing equation ( 22) by iv 1 + jv 2 on both sides and replacing f (i, j; t) by (iv 1 + jv 2 )C(i, j; t), we get the discrete population balance equation ( 9) after simplification.This shows that our Monte Carlo method is consistent with the population balance equation.4. Results and discussion.

PB model validation.
To validate the preliminary model developed in this paper, with both ε and φ values derived from untreated melanoma and PMN cases, we analyzed both the PMN melanoma aggregation results from parallel plate experiments and results from the model simulations in which the PMNs were treated with blocking antibodies.The average of the φ values calculated from the controlled (untreated cell cases) experiments was used in the coagulation kernel (equation 8), and the percentage of aggregates was recalculated for antibody blocked cases (Table 2).The results (Table 2) show that the simulations predict the percent aggregation within the experimental range in nearly all the cases.The experimental range was calculated as the mean minus standard deviation (SD) to the mean plus SD from three separate experiments.In one case (Anti-CD11a with shear rate 200/s), the range is not reported because the experimental data was the same for the three different experiments; therefore, no standard deviation could be calculated.If we take a range using comparable standard deviation as in the other cases, the simulation result in this case would be consistent to such a range.Thus, the overall consistency of the simulation results with experimental data shows the feasibility of using population balance models to simulate PMN melanoma aggregation in the near wall region of the parallel plate chamber.
While the numerical results are encouraging, it also demonstrates that to improve the accuracy of the model, a simple normalization of the initial population is not enough.For example, for the case of Anti-CD11a with shear rate 100/s, the simulation result is slightly out of the experimental range, in addition to the mismatch for the case with 200/s shear rate if we do not borrow a similar range of data In the current population balance model, we used the scaling invariance to introduce a normalization factor to suitably adjust the local concentration in the near wall region since it can not be measured accurately experimentally.However, the limitation of this technique is that we are not able to change the ratio of PMN and tumor cell concentrations.For example, suppose the ratio of PMN and melanoma cell concentration is 1 : 1 in the inlet; it may be changed to totally different ratio in the near wall region.
In order to estimate the cell concentrations in the near wall region directly, careful CFD simulations of cell transport within the chamber are under development [13].Such simulations can provide more precise quantification of the local cell number densities.Therefore our future approach is to perform direct numerical simulation of a statistically significant number of cells transported from a uniform reservoir concentration.Some preliminary study provides feasibility of employing CFD to perform a detailed parametric study across a range of melanoma and PMN reservoir concentrations, Reynolds numbers and Capillary numbers to obtain statistically relevant local concentration predictions for integrating with PB modeling.In our future works, such CFD simulations will be combined with the population balance modeling to provide more accurate determination of the local cell concentration.Such studies will also help our understanding of the more physiologically relevant case where the number of PMNs is significantly greater than melanoma cells (normal physiological levels are 2x10 9 − 7.5x10 9 PMNs per liter of human blood).In particular, it would be important to develop a PB model in the case where the ratio of PMNs to melanoma cells is high in the free stream.It is then necessary to take into account changes in spatial ratios of PMNs and melanoma cells populations from the free stream ratios to final surface ratios near the EC substrate, by using computational fluid dynamics simulations.

Figure 1 .
Figure 1.Cell aggregation in the near wall region.

3. 3 .
Determining the coagulation kernel function.The coagulation kernel presented in equation (3) takes the presence of cells into account through the introduction of the adhesion efficiency ε.Still, this formulation assumes a uniform shear field and a spatially homogeneous distribution of cells considered as rigid particles.

Figure 2 .
Figure 2. The average shear rate in the near wall region.

4. 3 .
Acknowledgments.This work was supported in part by the National Science Foundation grants NSF-DMS 0409297, NSF-DMR ITR-0205232 and NSF-CCF 0430349 (Q.D.), the National Institutes of Health grant CA-97306 (C.D.) and the National Science Foundation grant BES-0138474 (C.D.).The authors appreciate the support of the Penn State General Clinical Research Center (GCRC) funded by the NIH Grant M01-RR-10732.

Table 1 .
Adhesion efficiency values for P +P and T +P collisions for three shear rates.T + T was approximately zero.
3.5.Determining the normalization factor.The normalization factor was determined by first solving the population balance equation (

Table 2 .
Comparison of simulated tumor PMN aggregation results to experimental results that used PMNs treated with anti-CD11a or anti-CD11b.One possibility is that more realistic PMN melanoma ratios near the EC substrate are needed, in addition to taking into account the change of total population near the substrate.4.2.Discussions and future work.