Numerical simulation of the viral entry into a cell driven by receptor diffusion
NNumerical simulation of the viral entry into a celldriven by receptor diffusion
T. Wiegold , ∗ , Sandra Klinge , R. P. Gilbert , G. A. Holzapfel , Institute of Mechanics, TU Dortmund University,Leonhard-Euler-Strasse 5,44227 Dortmund, Germany Chair of Structural Mechanics and Analysis, TU Berlin, Strasse des 17. Juni135, 10623 Berlin, Germany Department of Mathematical Sciences, University of Delaware, 312 EwingHall, Newark, Delaware 19716, USA Institute of Biomechanics, Graz University of Technology, Stremayrgasse16/2, 8010 Graz, Austria Department of Structural Engineering, Norwegian University of Science andTechnology (NTNU), 7491 Trondheim, NorwayThe formal publication of this article can be found via:https://doi.org/10.1016/j.camwa.2020.12.012
Abstract
The present study focuses on the receptor driven endocytosis typicalof viral entry into a cell. A locally increased density of receptors at thetime of contact between the cell and the virus is necessary in this case.The virus is considered as a substrate with fixed receptors on its surface,whereas the receptors of the host cell are free to move over its membrane,allowing a local change in their concentration. In the contact zone themembrane inflects and forms an envelope around the virus. The createdvesicle imports its cargo into the cell. This paper assumes the diffusionequation accompanied by boundary conditions requiring the conservationof binders to describe the process. Moreover, it introduces a conditiondefining the energy balance at the front of the adhesion zone. The latteryields the upper limit for the size of virus which can be engulfed by thecell membrane. The described moving boundary problem in terms of thebinder density and the velocity of the adhesion front is well posed andnumerically solved by using the finite difference method. The illustrativeexamples have been chosen to show the influence of the process parameterson the initiation and the duration of the process.
Keywords : cell mechanics, receptor diffusion, moving boundary problem, finitedifference method ∗ Corresponding author, email address: [email protected] a r X i v : . [ q - b i o . CB ] J a n Introduction
The intense study of cell mechanisms has provided an important insight intothe uptake of various substances into a cell including viruses. Diagnosis andtherapy of diseases has reached a state in which nanomedicine concerned withdevices of nanoscale size is applied. With these, they can deliver low molecu-lar mass compounds, proteins and recombinant DNAs to focal areas of disease.Some examples are polymeric micelles, quantum dots, liposomes, polymer-drugconjugates, dendrimers, biodegradable nanoparticles, silica nanoparticles, etc.researched in laboratories, undergoing preclinical development, or already usedin hospitals [1, 2]. In addition, numerical methods in biomechanics and biomath-ematics become more relevant. A variety of numerical methods ranging frommultiscale finite elemente [3] and isogeometric shell formulations [4] to relaxationdynamics [5] find their application in biomechanics.The cell is surrounded by a plasma membrane which acts as the interfacebetween the cell and its surrounding environment. However, the membrane isnot absolutely impermeable and a transport of particles through the membraneis still possible. Among many different mechanisms, the most common processfor this purpose is the so-called endocytosis [6]. The main focus of the investi-gation of endocytosis has so far been on clathrin-mediated endocytosis (CME).During the CME, proteins create clathrin-coated pits which eventually buildwhole vesicles [7, 8].Various aspects of the endocytosis process are investigated in different chem-ical and biochemical contexts. Amongst others, the total internal reflection illu-mination with fluorescence correlation spectroscopy is used to measure ligand-receptor kinetic dissociation rate constants [9]. Trafficking phenomena are stud-ied based on internalization experiments utilizing multiscreen assay systems [10].Furthermore, biomimetic systems of lipid vesicles or supported bilayers with avariety of binder molecules deal with soft adhesion mediated by mobile bindersas shown in [11, 12]. The behavior and stability of adhesion complexes arealso addressed by using cell doublet [13, 14] and vesicles adhered to supportedbilayers [15].The experimental progress has parallely led to multiple theoretical models.A mathematical framework, based on probabilities for binding rates, introducinga random and a sequential driving mechanism for the receptors is provided in[16]. Discrete stochastic models for specific receptor-ligand adhesion as well asnon-equilibrium continuum models for the competition between different modesof junction remodeling under force are developed in [17], whereas a statisticalthermodynamic model of viral budding is presented in [18]. Several analyticalmodels for the endocytosis process utilize the description of the Stefan problem,[19, 20] and propose a solution relying on the error and complementary errorfunctions. These methods are particularly applied to the HIV-Virus [21] andSemliki Forest virus [20] but also for the uptake of nanoparticles [22]. Alterna-tively, Tseng and Huang [23] use the immersed boundary method to simulatethe endocytosis and to investigate the resistance of the water film in the contactarea. 2he current contribution uses the model presented in [19, 20] as a basis,however, it proposes an alternative form of the Stefan supplementary conditionwhere the focus is only set on the energetic aspects of the front itself, and aconsideration of the dissipation associated with the receptor transport alongthe cell membrane is not needed. Different from the previous models, the newformulation yields the upper limit for the size of virus able to enter the cell. An-other focus of the paper is the purely numerical solution of the moving boundaryproblem where no additional assumptions typical of analytical solutions such asthe speed factor [19, 20] are necessary. The chosen interpretation of the problemis advantageous for the fast simulation of various scenarios regarding the processparameters and their influence on the initiation and duration of the process. Asa final objective, the paper also introduces the notion of cooperativity [24, 25]into the numerical model and performs a more sophisticated study of the effec-tive receptor distribution on a virus and of its influence on the viral entry intoa cell.This contribution is structured as follows: A general overview of the uptakeprocess is summarized in Sect. 2, whereas Sect. 3 recapitulates the main aspectsof the free energy characteristic for the endocytosis process. After this intro-ductory part, Sect. 4 focuses on the definition of boundary and supplementaryconditions accompanying the driving diffusion differential equation. The bound-ary conditions define the flux balance, and supplementary condition representsthe energy balance at the adhesion front. Subsequently, Sect. 5 discusses thenumerical implementation for the 1D case and helical viruses, which is followedby a non-dimensional description of the problem in Sect. 6. The paper alsodeals with a common case of a spherical virus (Sect. 7) and extends the basicmodel by introducing the notion of cooperativity (Sect. 8). Finally, Sect. 9gives a comprehensive overview of results and provides a comparison with theexperimental observations and numerical results from the literature. The paperfinishes with a conclusion and an outlook.
In a basic view of mechanical adhesive contact between elastic surfaces, twophenomena which have a considerable influence on the underlying process coun-teract each other. A reduction in the free energy when surfaces with bondingpotential come into contact benefits the process, whereas an increase in freeenergy due to elastic deformation required to fit their shapes counteract theprocess. In the classic Hertzian theory of elastic deformation [26], two bodiescoming into contact deform in the contact area in such a way that they perfectlyfit. According to this approach any surface interactions such as Van der Waalsforces, which are induced by charge polarization in electrically neutral moleculesin close proximity, are excluded. However, these ’non-material’ effects have asignificant influence on the direct contact interaction [27]. This is illustrated byconsidering small elastic objects consisting of crystalline materials processed ina controlled environment. Such crystals show the appearance of unfulfilled or3angling chemical bonds distributed over a free surface. Bringing such objectsinto contact reduces the free energy of the system by forming bonds between thetwo surfaces. The objects joined in this way will not separate without additionalwork. Hence, not only compressive traction due to bulk elasticity but also an ad-hesive or tensile traction contributes to the contact [19]. General investigationsof the mechanisms on adhesive contact are presented in [28, 29, 30, 31].The same effects, both attractive and resisting interactions, appear in theadhesive contact of biological cells. However, due to their characteristic prop-erties compared to engineering materials, significant differences occur in thiscase. Having a remarkably lower elastic modulus than engineering materialsweakens the influence of the effect of elastic energy variations during contact[19]. Furthermore, cells are characterized by a fluid-like in-plane behavior[32].This enables the receptors of the cell to move within its membrane, enablingnew methods of incorporating free energy variations in the modeling of adhesivecontact.For a long time, electron microscopy has been used to provide valuable in-sight into the architecture of viruses [33]. Furthermore, fluorescent-labelingof viruses and cellular structures combined with fluorescence microscopy yieldmore dynamic information for the tracking of a single virus in live cells [34].In order to obtain 3D geometrical and distributional information, electron to-mography has also shown to be a powerful tool [35]. Thus, a large amount ofinformation on the architecture of viruses is already available. For example, aspike protein density of approximately 2800 µ m − is identified for Sars-Cov [36]and the value of approximately 5200 µ m − for an Alphavirus [18]. With regardto the geometry, the investigations have shown that the helical and sphericalvirus forms are predominant [37].In order to depict the process of viral entry into a cell, the situation presentedin Fig. 1 is considered. This 1D situation is suitable to simulate the endocytosisof helical virus into a cell. However, an extension to the rotationally symmetriccase and to the simulation of a spherical virus is straightforward (Sect. 7). In theinitial state, the virus has not yet reached the cell surface (Fig. 1a). Upon firstcontact, the virus gradually connects to the cell (Fig. 1b). In order to establisha connection between the virus and the cell, a generic repulsion between theirsurfaces needs to be overcome.The connection by binding receptors of the virus to receptors of the cellreduces the internal energy of the system. Upon completing a single receptor-ligand bond, the internal energy is reduced by k T C b , where k is the Boltzmannconstant, T the absolute temperature and C b the binding energy coefficient.The quantity driving the uptake process is receptor density ξ . Initially, thedensity of receptors on the cell surface amount to ξ and the correspondingcounterpart, the receptor density on the virus surface, amounts to ξ rv . In gen-eral, it holds that the density of receptors on the virus is larger than the one onthe cell surface and that the virus receptors are fixed, whereas the receptors ofthe cell are free to move across the membrane. Upon contact, receptors of thecell diffuse over the surface, connect to the receptors of the virus and build anenvelope around the virus. At the end of the process the envelope is closed over4 a) (b) Figure 1: (a) Initial configuration of the cell surface and the virus in a 2D setup.(b) State during the uptake where the virus is partially connected to the cell.the virus which has fully entered the cell.As opposed to metallic or covalent bonds, the bonds created during biologicaladhesion are relatively weak. Since the cell receptor density is the lower one, ingeneral, it dictates the amount of reduction in the internal energy. Typically,the resisting potential due to generic repulsion exceeds the reduction in internalenergy of the initial configuration of the system for a unit area of the membraneat ξ . Therefore, additional influences facilitate the creation or dissolving ofchemical bonds. Possible influences are catalytic agents, small temperaturechanges and small mechanical forces. It appears that a local change in receptordensity is necessary in order to create an adhesion zone between the virus andthe cell. An increasing local receptor density results in a greater reduction in thefree energy by completion of each additional bond. When the cell receptors andthe virus receptors are close to each other, a permanent interaction is presentdue to thermal stimulation.In the framework of chemical rate theory two distinct cases are differentiated[38]. In the area where ξ < ξ eq holds, the rate of bond breaking exceeds the rateof bond forming, so that no adhesive contact is established. In the area where ξ > ξ eq holds, the rate of bond forming exceeds the rate of bond breaking, anda strong adhesive contact is established. Condition ξ = ξ eq defines the chemicalequilibrium of the bonding reaction as well as the lower limit for the adhesionto start. Whereas it is known that ξ eq > ξ in the case of cell-virus contact,different assumptions can be made for the exact value of ξ eq . In a limitingcase, the chemical equilibrium requires all receptors of the virus to be bondedto the cell membrane such that ξ eq = ξ rv . It is more likely to expect that theequilibrium density ξ eq is in the range [ ξ , ξ rv ] and that it might vary during theprocess. Some of these aspects are considered in Sect. 8 on cooperativity, aneffect significantly influencing the exact value of ξ eq . However, the initial studyis performed for the most restrictive case ξ eq = ξ rv , which does not influencethe generality of the model.As an illustration, a schematic distribution of receptors over the cell and thevirus is depicted in Fig. 2a, whereas the corresponding density profile is shownin Fig. 2b. In the adhesive zone, the receptor density is constant and amountsto ξ eq . The density starts to grow outside the adhesive zone and, only when5 a) (b) ξ (cid:2) µ m − (cid:3) a ( t ) ξ eq ξ ξ + x [ µ m] Figure 2: (a) Schematic diagram for the receptor distribution over the cell andvirus surface. (b) Typical density profile.far away from it, moves towards the initial density value of the cell ξ . Sincethe size of the cell is magnitudes larger than the virus, we assume the receptordensity far away from the adhesion front to stay constant lim x →∞ ξ ( x, t ) = ξ .Consequently, the flux j of the receptors vanishes so that lim x →∞ j ( x, t ) = 0.One more peculiarity of diagram 2b is the front of the adhesion zone wherea jump of the receptor density occurs. The position of the front is denoted bya time dependent function a ( t ). The discontinuous profile is expected since theequilibrium density is high and certainly larger than ξ , whereas the receptordensity ahead of the front must be lower than ξ in order to stipulate the receptordiffusion. The values typical of the adhesion front play an essential role in themodel to be described, and are denoted by the subscript + in the subsequenttext. For example, ξ + denotes the receptor density at the front.The previous explanation shows that the whole process is regulated by thediffusion of receptors over the cell surface and their gathering in the adhesionzone. Accordingly, the motion of the receptors will be described by the diffusiondifferential equation ∂ξ∂t = − ∂j∂x (1)which states that the change of receptor density in time has to be equal to thenegative spatial change in the flux. Furthermore, following Fick’s first law, thereceptor flux j is proportional to the gradient of density, i.e. j = − m ∂ξ∂x , (2)such that the implementation of Eq. (2) into Eq. (1) leads to the alternativeexpression of the diffusion equation ∂ξ∂t = m ∂ ξ∂x . (3)6his equation defines the relation between the temporal and the spatial changesof the receptor density weighted by the mobility parameter m . Its evaluationgives insight into the evolution of receptor density for every point in front of theadhesion zone a ( t ) < x < ∞ . Equation (3) is a partial differential equation ofsecond order and requires additional boundary conditions in order to determinethe complete particular solution. These two conditions will be defined in theupcoming sections. In order to define the free energy characteristic of the simulated process, thesystem including a large number of receptors is treated analogously to the caseof an ideal gas with a large number of non-interacting particles N . In such acase, the entropy of a single particle, belonging to a system in equilibrium, isexpressed by k ln[( A/ Λ ) ( c/N )] [19]. Here, A is the considered surface, N/A is the areal density ξ , c is a numerical factor and Λ a molecule length scale[39]. However, the latter two quantities ( c and Λ) do not play any role for thedescription of our process since it does not depend on the absolute entropy buton its change. This change is described with respect to the initial state of the cellwith uniform density ξ which is chosen to be the reference state. The relativeentropy of a single receptor at density ξ is then described by its difference tothe reference state and is calculated according to [19] by k ln (cid:18) cξ Λ (cid:19) − k ln (cid:18) cξ Λ (cid:19) = − k ln (cid:18) ξξ (cid:19) . (4)With the expression for one receptor at hand, and assuming that the gradientin the distribution is small and that the local distribution is indistinguishablefrom an equilibrium distribution at local density, the free energy E e per unitarea of membrane surface associated to the receptor distribution at absolutetemperature T turns into E e = k T ξ ln (cid:18) ξξ (cid:19) . (5)Moreover, the chemical potential χ is defined as the local change in the freeenergy per receptor, χ ( x, t ) = ∂E e ∂ξ = k T (cid:20) ln (cid:18) ξξ (cid:19) + 1 (cid:21) . (6)Finally, the mean receptor speed is assumed to be proportional to the spatialgradient of the chemical potential, i.e. v r = − mk T ∂χ∂x = − mξ ∂ξ∂x , (7)where the motion of the receptors is controlled by the mobility parameter m .7 Boundary and supplementary conditions
The full description of the adhesion front motion relies on a problem formula-tion including the differential equation (3) along with two boundary conditionsdescribing the flux at the ends of the domain and along with a supplementarycondition.
The boundary conditions on the unbound area are concerned with the quanti-tative description of the flux of receptors through the adhesive front. Followingthe Leibniz integration rule of the global form, this condition is derived fromEq. (1) as ( ξ eq − ξ + ) v + + j + = 0 , (8)or by using Fick’s first law as( ξ eq − ξ + ) v + − m (cid:20) ∂ξ∂x (cid:21) + = 0 . (9)Here, the first term denotes the amount of receptors required for the advance-ment of the front, and the second term denotes the amount of receptors providedby the flux. Both previous formulations depend on the front velocity definedin terms of the front position v + = d a d t . Equation (9) is consistent with theassumption (7), which can easily be shown as follows. First, the flux is assumedto be proportional to the receptor distribution ξ and the mean receptor velocity v r : j = ξ v r . (10)By incorporating Eq. (7) into Eq. (10), the flux turns into j = − mk T ξ ∂χ∂x = − m ∂ξ∂x , (11)as predicted by Fick’s first law. The supplementary condition is provided by considering the energetic aspectsof the front motion. The change of the receptor distribution as well as of themembrane shape leads to several contributions to the free energy of the system.However, the crucial observation is that the difference in the energy ahead ofand behind the front results in the front motion, which is expressed as follows, E − − E + = E kin . (12)Here, E − denotes the energy behind the adhesion front, E + is the energy aheadof the front and E kin is the kinetic energy of the front itself.8he term related to the energy behind the front is built of three contribu-tions, all denoted by superscript − , E − = E − b + E − e + E − κ . (13)These terms have the following physical meaning: E − b is the energy relatedto the binding of receptors, E − e the energy related to the entropy and E − κ theenergy related to the bending of the membrane. The reduction in the free energydue to the binding of receptors of the cell to receptors of the virus is defined asfollows E − b = − k T C b ξ eq . (14)This term is proportional to the reduction of energy caused by a single bondbetween two receptors − k T C b , and to the total amount of created bonds ξ eq dictated by the virus. As stated in [19], the binding energy coefficient C b typ-ically takes values in the range 5 < C b <
35. The second term describes theenergy associated with the entropy of receptors E − e = k T ξ eq ln (cid:18) ξ eq ξ (cid:19) , (15)which is required to bring the density from its reference value ξ to the densityof the virus ξ eq . This term will result in an increase in the free energy since itholds ξ < ξ eq . The third term of (13) is concerned with the bending of themembrane caused by the geometry of the virus E − κ = 12 k T B κ . (16)Here, a simplified case is considered corresponding to the theory of the bendingof a plate. Factor B represents the non-dimensional numerical parameter forthe bending stiffness, which is in the range of 10 to 30 and κ = 1 /R representsthe curvature, which is constant for a spherical virus and which depends on theradius of the virus R . Thus, the whole energy behind the front is then definedby the expression E − = − k T C b ξ eq + k T ξ eq ln (cid:18) ξ eq ξ (cid:19) + 12 k T B κ . (17)In the second step we consider the energy ahead of the front, denoted by su-perscript +. Binding between the cell and the virus exclusively takes place inthe area behind the front and thus does not have any influence on the energyahead of the front. However, corresponding parts E +e , the energy related to theentropy and E − κ , the energy related to the curvature of the membrane remainavailable. Moreover, a term E +v , the energy related to the motion of receptorsalso has to be taken into consideration. In summary, the following terms canbe counted ahead of the front: E + = E +e + E + κ + E +v . (18)9n the present contribution, we assume that the curvature behind the front ismuch smaller than the one caused by the contact with the virus. This justifiesthe assumption of a vanishing influence to the energy associated to the bendingof the membrane E + κ = 0 . (19)The energy ahead of the front related to the entropy is expressed in the sameway as the energy behind the front as E +e = k T ξ + ln (cid:18) ξ + ξ (cid:19) . (20)It describes the energy needed in order to bring the initial receptor density ξ tothe value ξ + . Contrary to the contribution behind the front, this term resultsin a reduction of the free energy since ξ > ξ + . The contribution due to themotion of the receptors reads E +v = 12 m r ξ + v = 12 m r m ξ + (cid:18) ∂ξ + ∂x (cid:19) . (21)It represents the kinetic energy of all receptors ahead of the front moving towardsthe front with their corresponding velocity v r and the cell receptor mass m r .With Eqs. (19) - (21), the total energy ahead of the front is defined by E + = k T ξ + ln (cid:18) ξ + ξ (cid:19) + 12 m r m ξ + (cid:18) ∂ξ + ∂x (cid:19) . (22)Finally, the difference between the energies of the two sides of the front actsas driving force for the front motion. The kinetic energy of the front is thencharacterized by the mass of the front m rr ξ eq and the front velocity v + E kin = 12 m rr ξ eq v . (23)Here, m rr represents the mass of a receptor pair including the cell receptor andthe virus receptor which are bonded to each other. Combining Eq. (17), (22)and (23) leads to the expression − ξ eq C b + ξ eq ln (cid:18) ξ eq ξ (cid:19) + 12 Bκ − (cid:34) ξ + ln (cid:18) ξ + ξ (cid:19) + 12 m r k T m ξ + (cid:18) ∂ξ + ∂x (cid:19) (cid:35) = 12 ξ eq m rr k T v , (24)which is the final form of the supplementary condition and closes the formulationof the moving boundary problem. In summary, the change of the receptor distribution is described by a systemof differential equations consisting of (3), (9) and (24). The finite difference10ethod has been chosen for the solution of the underlying system of differentialequations. According to this approach, all derivatives are replaced by expres-sions dependent on discrete values of the function for the nodes of a chosenlattice. Thus the differential equations are transformed into a system of alge-braic equations. An implicit scheme is used, with the following approximationsfor the derivatives ∂ξ∂t = ξ j +1 i − ξ ji ∆ t , ∂ξ∂x = ξ j +1 i − ξ j +1 i +1 ∆ x , ∂ ξ∂x = ξ j +1 i − − ξ j +1 i + ξ j +1 i +1 ∆ x . (25)Here, subscript i denotes the spatial position and superscript j denotes the time.The implementation of relationships (25) into the system (3), (9) and (24) leadsto the following discretized formulation of the problem: ξ j +1 i − ξ ji ∆ t = m ξ j +1 i − − ξ j +1 i + ξ j +1 i +1 ∆ x , i = 1 , ..., p, j = 1 , ..., n, (26a) (cid:104) ξ eq − ξ j +1+ (cid:105) v j +1+ + m ξ j +1+ − ξ j +11 ∆ x = 0 , (26b) (cid:2) E − (cid:3) − ξ j +1+ ln (cid:32) ξ j +1+ ξ (cid:33) + 12 m r k T m ξ j +1+ (cid:32) ξ j +11 − ξ j +1+ ∆ x (cid:33) − (cid:20) ξ eq m rr k T v j +1+ 2 (cid:21) = 0 . (26c)In Eq. (26a), variable p refers to the total number of points ahead of the front,except the last point where the influence of the flux vanishes and where thereceptor density is kept at the initial value ξ . Variable n refers to the totalnumber of time steps. Furthermore, the conditions given in Eqs. (26b)-(26c)are valid at the front. In Eq. (26c), E − is an abbreviation for the contributiondefined in (17). This term does not depend on the density ξ i and quantities atthe front ξ + and v + , and thus represents a constant during the process.The solution of system (26) yields values for receptor densities ξ j +1+ and ξ j +1 i , i = 1 , ..., p, and front velocity v j +1+ . The latter is, in postprocessing, usedto evaluate the front position in the incremental form: a j +1 = a j + v j +1+ ∆ t ,where j + 1 and j are two subsequent time steps. The numerical examples chosen simulate the process of virus uptake into the cell.In the simulations, it is assumed that a helical virus of size D = 0 . µ m, comesinto contact with a much larger cell such that the cell curvature is negligible(Fig. 1). Due to the axial symmetry of the virus, the problem is treatedin a 1D representation which assumes the unit width of the active domain.The initial density of cell receptors is set to ξ = 1000 µ m − , whereas theinitial density of virus receptors is set to ξ eq = 4800 µ m − . Time increment∆ t = 1 e − s and space increment ∆ x = 1 e − µ m are used for the numericalsimulations. An overview of the chosen process parameters is given in Table 1.11 aterial parameters Receptor density on cell surface ξ µ m − Receptor density on virus ξ eq µ m − Receptor mass m r
400 kDaMass of a receptor pair m rr
800 kDaBinding energy coefficient C b − Numerical bending stiffness parameter B − Curvature of the virus κ µ m − Mobility parameter m . µ m / sVirus diameter D . µ mTable 1: Process parameters used in simulations.These values belong to the corresponding admissible ranges and are also usedin [19]. The convergence of results has been checked by varying the time andthe space increment. Time increment ∆ t has been decreased in the interval∆ t = 1 e − − e − s which has caused a change of results in ξ + for maximally0 . x in range 1 e − − e − µ m causedthe changes in ξ + up to 2 . m =1 µ m /s. The density profiles for different time steps during the simulation arepresented in Fig. 3. The diagrams show a fast decrease in receptor density,particularely at the beginning of the process. After 300 time steps, the densityat the front only amounts to ≈
50% of its initial value. This rapid decline indensity at the front slows down in the course of the further process. R ece p t o r d e n s i t y ξ [ µ m − ] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m]Time step 240 Time step 480 Time step 720 Time step 860 Time step 1200 Figure 3: Receptor density ξ over the cell surface x for the first 1200 time steps,corresponding to a simulated time of 0 .
12 s.Figure 4 monitors the advancement of the front and the position of the virusduring its entry into the cell in 1D representation. The position of the virus12s related to the position of the front through length a , determining the sizeof contact area. For a helical virus, a 3D visualization is also possible due tothe axial symmetry, as shown in Fig. 5 which compares the endocytosis of avirus into cells with different receptor mobilities. In the top row, the mobilityis set to m = 1 µ m /s, whereas half of this value m = 0 . µ m /s is used insimulations in the lower row. Naturally, the first process is faster and the viralentry is accomplished earlier than in the second case. The increasing numberof time steps between the three states indicates the gradual decrease and finalstagnation of the velocity of the process, an issue also studied in the followingexample. R ece p t o r d e n s i t y ξ [ µ m − ] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] a a = πR a = πR a = πR a = πR Time step 146 Time step 385 Time step 743 Time step 1228
Figure 4: Visualization of the front motion and of the formation of the envelopearound the virus with diameter D = 0 . µ m.The governing equation (3) of the process depends on a single process pa-rameter, namely on mobility m . The parameter represents a measure for thecapability of receptors to move over the cell surface, and thus is in a direct cor-relation with the amount of receptors provided for the adhesion with the virus.The influence of the mobility on the velocity of the front and on the receptordensity has been studied on the basis of a set of simulations, as shown in Fig. 6.Here, mobility parameter m has been varied in the range [0 . µ m / s-1 µ m / s].Figure 6a shows the dependence of velocity v + on the mobility and clearly con-firms the rapid decrease in velocity at the beginning followed by a stagnation,as already observed in the previous test (Fig. 4). The value of the mobilitydoes not affect the form of the velocity diagrams. However, a higher velocitycorresponds to a higher mobility. This observation is in agreement with thephysical character of the mobility describing the ability of receptors to movetowards the adhesion zone. For lower values of m , fewer receptors are providedto connect the cell with the virus. Therefore, the evolution of the adhesion zoneand the velocity of the front are slowed down. An analogous trend is observed13 = 214 e − s a = 0 . µ m t = 609 e − s a = 0 . µ m t = 1227 e − s a = 0 . µ m a = 0 . µ m a = 0 . µ m a = 0 . µ m Figure 5: Endocytosis of a helical virus. Comparison of the process state at thesame time for two different mobilities. The upper row shows results correspond-ing to a higher mobility m = 1 µ m / s. Results presented in the bottom row arerelated to lower mobility m = 0 . µ m / s.for the dependency of the receptor density at the front on the mobility shownin Fig. 6b. . . . . . . F r o n t v e l o c i t y v + [ µ m / s ] Time [1 e − s] m = 1 . µ m / s m = 0 . µ m / s m = 0 . µ m / s(a) (b) 10008006004002000 0 400 8000 1200Time [ 1 e − s] F r o n t d e n s i t y ξ + [ µ m − ] m = 1 . µ m / s m = 0 . µ m / s m = 0 . µ m / s Figure 6: (a) Velocity of the adhesion front vs. time. (b) Evolution ofthe receptor density at the adhesion front. Mobility is varied in the range0 . µ m / s − µ m / s.An important influence on the process is also imposed by the fixed receptor14 . . . . . . . .
00 100 200 300 400 500 ξ eq = ξ eq = ξ eq =Time [1 e − s] F r o n t v e l o c i t y v + [ µ m / s ] (a) 3200 µ m − µ m − µ m − . . . . . . . . ξ = ξ = ξ =Time [1 e − s] F r o n t v e l o c i t y v + [ µ m / s ] (b) 1500 µ m − µ m − µ m − Figure 7: (a) Evolution of the velocity of the adhesion front for different densities ξ eq . Density ξ is set to 1000 µ m − . (b) Evolution of the velocity at the adhesionfront front for different initial densities ξ . Density ξ eq is set to 4800 µ m − .density ξ eq of the virus, initially chosen to dictate the amount of receptorsrequired for the virus-cell connection. The velocity of the adhesion front v + fordifferent values of ξ eq is shown in Fig. 7a. Here, the receptor density of the cellis set to ξ = 1000 µ m − and the mobility is set to m = 1 µ m /s. The formof the velocity diagrams does not change, although the different constellationsare taken into consideration. The velocity of the adhesion front v + decreaseswith increasing density ξ eq , which is to be expected since a larger number ofreceptors is necessary in order to achieve a front advancement.Similar simulations are conducted for different values of the initial receptordensity ξ , while the receptor density of the virus is set to ξ eq = 4800 µ m − (Fig. 7b). Again, the initial configuration does not affect the form of the dia-grams, whereas a larger density ξ corresponds to higher velocities. The requiredamount of connected receptors has been fixed at a constant value in all the sim-ulations. However, only the initiation of the process requires a higher amountof bonds. Once contact between the cell and the virus has been established, thenumber of necessary receptors decreases. The amount of bonds required for thecontact between the cell and the virus cannot fall below a minimum value. Theevolution of the required cell receptor density can be easily implemented in thedeveloped code by assuming ξ eq to be a function of time. The simulations inthis case (results not shown here) indicate an accelerated viral entry into thecell as a consequence of the decrease of the required receptor density.15 Non-dimensionalization
In some physical systems, non-dimensionalization is applied to suggest thatit is more convenient to measure certain quantities relative to an appropri-ate unit. These units refer to quantities intrinsic to the system. The non-dimensionalization procedure relies on replacing dimensional quantities by thenon-dimensional ones within the differential equation and within the correspond-ing boundary and supplementary conditions. An important advantage of thenon-dimensional analysis is that it reduces the number of relevant process pa-rameters and thus facilitates the parameter study as performed in the previoussection.The present model for the viral entry includes following dimensional quan-tities which can be expressed in terms of their non-dimensional counterparts x = l s x (cid:48) , t = τ s t (cid:48) , ξ = ξ s ξ (cid:48) . (27)Here, l s , τ s , ξ s are properly chosen scaling parameters for space, time and density,and the prime symbol denotes the dimensionless quantities. The introduction oftransformation (27) in differential equation (3) yields its non-dimensional form ∂ξ (cid:48) ∂t (cid:48) = m τ s l ∂ ξ (cid:48) ∂x (cid:48) . (28)The same procedure can now be applied to the flux-boundary condition (9) andthe supplementary energy condition (24) which, amongst others, depend onfront velocity v + . The latter intrinsically includes the derivative with respectto time, which yields v + = ∂a∂t = ∂ ( l s a (cid:48) ) ∂ ( τ s t (cid:48) ) = l s τ s v (cid:48) + . (29)Bearing in mind transformations (27) and intermediate result (29), the fluxboundary condition and the energetic supplementary condition have the non-dimensional form ( ξ (cid:48) eq − ξ (cid:48) + ) v (cid:48) + − m τ s l (cid:20) ∂ξ (cid:48) ∂x (cid:48) (cid:21) + = 0 , (30) C (cid:48) − ξ (cid:48) + ln ξ (cid:48) + − m r k T m l ξ (cid:48) + (cid:18) ∂ξ (cid:48) + ∂x (cid:48) (cid:19) = 12 m rr k T l τ ξ (cid:48) eq v (cid:48) , (31)where the following abbreviations apply ξ (cid:48) + = ξ + ξ s , ξ (cid:48) eq = ξ eq ξ s , (32) C (cid:48) = − ξ (cid:48) eq C b + ξ (cid:48) eq ln ξ (cid:48) eq + 12 B κ ξ s . (33)16revious formulation (30)-(31) calls upon the introduction of additional non-dimensional process parameters m = m τ s l , m r = 12 m r k T m l , m rr = 12 m rr k T l τ ξ (cid:48) eq . (34)where m represents the dimensionless mobility, m r is the non-dimensional cellreceptor mass and m rr the non-dimensional mass of a receptor couple. Thisshort notation yields the final form of the non-dimensional moving boundaryproblem ∂ξ (cid:48) ∂t (cid:48) = m ∂ ξ (cid:48) ∂x (cid:48) , (35)( ξ (cid:48) eq − ξ (cid:48) + ) v (cid:48) + − m (cid:20) ∂ξ (cid:48) ∂x (cid:48) (cid:21) + = 0 , (36) C (cid:48) − ξ (cid:48) + ln ξ (cid:48) + − m r ξ (cid:48) + (cid:20) ∂ξ (cid:48) + ∂x (cid:48) (cid:21) = m rr v (cid:48) . (37) After deriving the non-dimensional formulation of the problem, a further im-portant step is choosing characteristic scaling parameters. Typically, these areadapted to the system properties. In the present case, half of the arclength ofthe virus is assumed as the characteristic length, namely l s = πR . On the otherhand, the characteristic time is chosen as the time necessary to complete thevirus uptake at a constant unit velocity v = 1 µ ms . Accordingly, the time scalingparameter reads τ s = l s /v . The initial cell receptor density is chosen as thelast scaling parameter, such that it holds ξ s = ξ .The results for a non-dimensional analysis are presented by examples inves-tigating the influence of non-dimensional parameters m, m r , m rr , C (cid:48) and ξ (cid:48) eq on the front velocity and front receptor density. To this end, first, the referencevalues for the scaling parameters are set as follows: l s = πR = 0 . µ m and τ s = l s /v = 0 . R = 0 . µ m is assumed as the referencevirus size. The density scaling parameter takes the value ξ s = ξ = 1000 µ m − .The results of the analysis parameter analysis are presented in Fig. 8. Theyshow that the dimensionless mobility m has an important influence on the frontvelocity (Fig. 8a), whereas the dimensionless mass m r mainly influences thefront receptor density (Fig. 8b). The variation of the dimensionless equilibriumdensity ξ (cid:48) eq has an important influence on both quantities (Fig. 8c). The effect of C (cid:48) is similar to the one of ξ (cid:48) eq , whereas the variation of the dimensionless couplemass m rr hardly affects the results (results not shown here). Amongst others,a higher mobility corresponds to the higher front velocity, whereas higher mass m r and equilibrium density ξ (cid:48) eq cause a higher front density.17 a) 2 . . . . . . . F r o n t v e l o c i t y v + / v [ − ] Time t/t [ − ]0 . . . . . m = 1 . m = 1 . m = 0 . . . . . . . . . F r o n t d e n s i t y ξ + / ξ [ − ] Time t/t [ − ]0 . . . . . m = 1 . m = 1 . m = 0 . . . . . . . . F r o n t v e l o c i t y v + / v [ − ] Time t/t [ − ]0 . . . . m r = 0 . m r = 0 . m r = 0 .
04 1 . . . . . . . F r o n t d e n s i t y ξ + / ξ [ − ] Time t/t [ − ]0 . . . . . m r = 0 . m r = 0 . m r = 0 . . . . . . . . F r o n t v e l o c i t y v + / v [ − ] Time t/t [ − ]1 . . . ξ ′ eq = 7 . ξ ′ eq = 4 . ξ ′ eq = 2 . . . . . . . . F r o n t d e n s i t y ξ + / ξ [ − ] Time t/t [ − ]1 . . . . ξ ′ eq = 7 . ξ ′ eq = 4 . ξ ′ eq = 2 . Figure 8: The non-dimensional analysis of the front velocity and the front den-sity. (a) The influence of non-dimensional mobility m . Remaining parametersare kept constant: m r = 0 . m rr = 0 . C (cid:48) = 7 . ξ (cid:48) eq = 4 .
8. (b) The influ-ence of non-dimensional receptor mass m r . Remaining parameters: m = 1 . m rr = 0 . C (cid:48) = 7 . ξ (cid:48) eq = 4 .
8. (c) The influence of non-dimensional equi-librium density ξ (cid:48) eq . Remaining parameters: m = 1 . m r = 0 . m rr = 0 . C (cid:48) = 7 .
53. 18
Rotationally symmetric case - Spherical virus
Whereas the previous analysis focuses on the 1D formulation applicable for thesimulation of helical viruses, the real situation is commonly concerned withthe rotationally symmetric geometry and spherical viruses. The problem for-mulation in this case is slightly different and requires an adaptation of thediffusion equation, whereas the boundary and supplementary conditions remainunchanged. The extension of the diffusion equation implies the introduction ofan additional term compensating for the radial dependency. Thus, the differen-tial equation turns into ∂ξ∂t − m ∂ ξ∂x − m x ∂ξ∂x = 0 , (38)where the last term on the left-hand side is the new contribution. The discretizedcounterpart of Eq. (38) is ξ j +1 i − ξ ji ∆ t − m ξ j +1 i − − ξ j +1 i + ξ j +1 i +1 ∆ x − m x ξ j +1 i +1 − ξ j +1 i − x = 0 , (39)where i = 1 , ..., p is the counter related to the spacial discretization and where j = 1 , ..., n is the counter corresponding to the time discretization. Variable x in the last term in (39) represents the distance from the first contact point andalso can be written in a discretized form as x = i ∆ x , which leads to a condenseddiscretization formulation ξ j +1 i − ξ ji ∆ t − m i ∆ x (cid:104) (2 i − ξ j +1 i − − iξ j +1 i + (2 i + 1) ξ j +1 i +1 (cid:105) = 0 . (40)Finally, the procedure explained in Sect. 6.1 provides the non-dimensinonalform of Eq. (38): ∂ξ (cid:48) ∂t (cid:48) − m ∂ ξ (cid:48) ∂x (cid:48) − m x (cid:48) ∂ξ (cid:48) ∂x (cid:48) = 0 . (41)The simulation of the virus uptake for a rotationally symmetric case is demon-strated by the example of the Alphavirus and by the process parameters sum-marized in Tab. 1. The achieved results are summarized in Fig. 9. They showthe 3D vesicle that is built during the process of the viral entry. Different tothe endocytosis of a helical virus, the process advances at a rather constantrate throughout the simulation. This change in behavior can be explained byproviding additional receptors due to the radial dependence. The same argu-mentation explains the higher velocity of the uptake process by a spherical viruscompared to the helical one. For the chosen parameters, the simulations predicta required time in the range for ultra-fast-endocytosis [40]. This fast behavioris expected, since the mobility has a rather high value. Viruses often connectto receptors with a lower mobility. By reducing mobility to m = 0 . µ m / s thetime increases and matches values for kiss-and-run-endocytosis [41]. The modelpredicts a shorter duration of the process than it is typical of the clathrin-mediated-endocytosis (15-20 s) [42]. 19 = , a = π R t = , a = π R t = , a = π R t = 1 , a = π R Figure 9: The 3D visualization of the endocytosis of a spherical virus withdiameter D = 0 . µ m. 20igure 10: Schematic representation of cooperativity during endocytosis (thefree receptors are not shown). Amongst others, cell adhesion deals with cooperativity, an effect which is ex-plained by considering a patch of unit length depicted in Fig. 10. As soon asreceptors create bonds, they smoothen out the surrounding membrane whichmakes it easier for additional receptors to create a bond and strengthens theadhesion between the virus and the cell membrane [24, 25]. This effect is knownas cooperativity. It has extensively been investigated experimentally and the-oretically. Different experiments are performed depending on the state of theadhesion process. The fluorescence recovery experiments are performed in or-der to analyze the equilibrated contact zone during the process, whereas themicropipette experiments are performed in order to analyze the initial contact.Lipid vesicles with anchored receptor molecules are often used in order to re-semble important aspects of cell adhesion.In order to study the binding cooperativity, two classes of numerical modelsare considered. The first class describes the membranes as continuous in spacewith continuous concentration profiles on the membrane [43]. The second classdescribes the membranes as discrete and the receptors as single molecules [44].Numerical solutions of the dynamic properties are studied by reaction-diffusionequations in the first class [45] of models and by Monte Carlo simulations in thesecond class [46]. The information obtained in such a way is complementary tothe model presented in this contribution.The cooperativity changes the amount of receptor bonds that will create anequilibrium state upon connection between the virus and the cell given in [24]according to ξ eq req = c κ b k T l K ξ ξ . (42)Here, ξ eq req is the required amount of receptors that need to bind in order tocreate adhesion between the virus and the cell. Symbol c denotes a dimension-less prefactor acquired from Monte Carlo simulations, usually ranging between10 - 15. Furthermore, the effective rigidity κ b can be calculated from the bend-ing rigidities of two apposing membranes as κ b = κ κ / ( κ + κ ). For thesimulations here, it is set to 40 k T . Quantity l we is the binding range depend-ing on the interaction range of the two binding sites, of the flexibility of theirmolecules and of the membrane anchoring. It describes the difference between21 a) (b) Small initial binding rangeHigh initial binding range Figure 11: (a) Contact with the virus with two kinds of receptors. (b) Spatialdistribution of different types of receptors on the virus membrane for two chosenconfigurations.the smallest and the largest local membrane separation at which the receptorscan bind. Quantity K pl is the two-dimensional equilibrium constant in the caseof two opposing planar, supported membranes within binding separation of thereceptor–ligand bonds.Two illustrative examples are performed in order to analyze the influenceof cooperativity for different binding ranges. As in the previous example theinitial receptor density of the cell is set to ξ = 1000 µ m − and the receptordensity of the virus is set to ξ eq = 4800 µ m − . The first group of simulationsconsiders a virus with its lower half covered by receptors with a smaller bindingrange and its upper half by receptors with a larger binding range. The lowerhalf is characterized by a binding range of l we = 1 nm resulting in the requiredreceptor density ξ eq req = 2265 µ m − , while the upper half is characterizedby a binding range of l we = 1 . ξ eq req = 3262 µ m − . An opposite situation is considered in the second groupof simulations presented in Fig. 11.Numerical results for the described examples are shown in Fig. 12. Thetransition between the areas with different receptor types manifests itself byeither a jump or a kink in the corresponding diagrams. The velocity is affectedmostly by the change of the required density. In the area with a smaller bindingrange less receptors are required, significantly increasing the velocity of theprocess. The diagrams for the second setup show similar results to the firstsetup, however, the change from the lower to the upper half is significantlydelayed. Here, the initial velocity is much higher in the first case such that thevirus is almost enclosed at the time step 400. Contrary to this, the velocity atthe end of the process is higher in the second case. Consequently, both virusesneed approximately 600 time steps for their entry into the cell. Exact valuesare 611 and 622 time steps for the first and second example respectively. Thevalues do not exactly match due to the different velocities at the beginning ofthe process and due to the transition between regions with different receptors.22 ime step 130 Time step 260 Time step 390 Time step 520 Time step 650 R ece p t o r d e n s i t y ξ [ µ m − ] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m](a) µ µ F r o n t v e l o c i t y v + [ µ m / s ] Time [1 e − s]Time [1 e − s]Time [1 e − s] F r o n t p o s i t i o n a [ µ m ] F r o n t d e n s i t y ξ + [ µ m − ] Time step 128 Time step 256 Time step 384 Time step 512 Time step 640 R ece p t o r d e n s i t y ξ [ µ m − ] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m] x -distance [ µ m](b) µ µ F r o n t v e l o c i t y v + [ µ m / s ] Time [1 e − s]Time [1 e − s]Time [1 e − s] F r o n t p o s i t i o n a [ µ m ] F r o n t d e n s i t y ξ + [ µ m − ] Figure 12: (a) Influence of the cooperativity for a smaller binding range in thelower half ξ eq req = 2265 µ m − and a larger one in the upper half ξ eq req = 3262 µ m − . (b) Influence of the cooperativity for a larger binding range in the lowerhalf ξ eq req = 3262 µ m − and a smaller one in the upper half ξ eq req = 2265 µ m − . Top row: Receptor density over the cell surface for different time steps.Bottom row: Velocity of the front, position of the front and receptor density atthe front over time. Chosen process parameters are K pl = 0 . e − and c = 13.23 . . . . . . N o r m a li ze d f r o n t p o s i t i o n a ( t ) / π / R [ − ] Normalized time t m/ ( π R ) [ − ](a) 0 . . . .
000 4 8 12 16 20 N o r m a li ze d f r o n t v e l o c i t y v ( t ) π R / m [ − ] Normalized time t m/ ( π R ) [ − ](b) Figure 13: Normalized front position (a) and normalized front velocity (b) versusscaled time.
The model developed gives insight into some specific features of the process andenables its profound analysis in the context of impeding and hindering the viralentry. Amongst others, it enables a study of the position of the front and itsvelocity during the process, an analysis of the admissible values for the radiusand of the duration of the process depending on different process parameters.The study presented uses the parameters listed in Tab. 1, if not stated otherwise.
The position of the front and its velocity are two characteristic indicators ofthe viral entry, giving insight into the current state of the process and enablingthe estimation of its total duration. The evolution of these indicators (Fig.13) shows that the front advances continuously during the process, whereas itsvelocity decreases with the strongest decline in the beginning, and an almostconstant value at the end of the simulation. Similar behavior is shown in thework by Freund und Lin [19]. However, while the overall trends in both worksare the same, some interesting phenomena can be identified with regard to thelimiting behavior. The current model experiences a weaker decline in the veloc-ity resulting in an almost constant velocity towards the end of the simulation,which yields a linear advancing of the front. However, the velocity in [19] movestowards zero which causes the process to experience almost no progress towardsthe end. 24 .2 Virus radius
Apart from the front position and velocity, the virus radius also gives importantinformation on the process, especially with regard to its initiation. This partof the analysis relies on the consequences of the energy balance (24). Since theright-hand side of this equation represents the kinetic energy, it directly followsthat the difference between the energy behind and ahead of the front on theleft-hand side has to be non-negative. In the beginning of the process, when thedensity distribution is uniform and the front has not yet been established, theenergy ahead of the front does not contribute to the total amount. Therefore,the part of the energy behind the front can be seen as an initial barrier that mustbe overcome in order to start the process. A study of the limiting case, wherethe front velocity approaches to zero, yields the expression for the maximalradius R max = (cid:114) B (cid:44) (cid:115) ξ eq C b − ξ eq ln (cid:18) ξ eq ξ (cid:19) (43)under the condition that C b − ln (cid:18) ξ eq ξ (cid:19) > . (44)By assuming the short notation for the receptor density ratio ˜ ξ = ξ ξ eq , the valuefor the critical density ratio results in˜ ξ crit = e − C b . (45)It is important to mention that the receptor density ratio is limited from bothsides. On one hand, it holds ˜ ξ = ξ ξ eq ≤ ξ ≤ ξ eq . On the other hand,it holds ˜ ξ > ˜ ξ crit due to the condition (44). The expressions (43) and (45)are now used to study the values of maximum radius. These results are shownin Fig. 14a, where the receptor density ratio is varied in the admissible semi-open range ( ˜ ξ crit , ξ crit according to Eq. (45). The correlation between the radius and the density ratiohas also been studied in the works by Gao et al. [20]. They derive similarexpressions for the limiting radius and the critical density ratio, however, theirmodel provides the expression for the minimal radius. In the present study aswell as in [20], the maximum radius increases rapidly as ˜ ξ approaches its criticalvalue. This quantity (43) also depends on the binding energy coefficient C b which is studied in Fig. 14b for three different values of the receptor densityratio. Here, the lower line corresponds to the higher ratio ( ˜ ξ = 1), and theupper line to the lower ratio ( ˜ ξ = 0 . C b , which becomes less pronounced asits value increases. A significant aspect of the virus entry is the duration of the complete processas well as its dependence on different process parameters. The viral uptake via25 a) 0 . . . . . . M a x i m u m v i r u s r a d i u s R m a x [ µ m ] Density ratio ˜ ξ [ − ]10 − − (b) 0 . . . . .
00 30 35 M a x i m u m v i r u s r a d i u s R m a x [ µ m ] Binding energy coefficient C b [ − ]5 10 15 20 25 ˜ ξ = 1 . ξ = 0 . ξ = 0 . Figure 14: (a) Maximum virus radius versus receptor density ratio ˜ ξ . Thevertical dashed line represents the critical value ˜ ξ crit . (b) The maximum virusradius versus the binding energy coefficient C b for different receptor densityratios ˜ ξ .endocytosis ranges through different time scales. Ultra-fast-endocytosis takes50-300 ms [40], while kiss-and-run-endocytosis takes approximately 1 s [41].In the presented model, several parameters have a significant influence on therequired time for the process. Three major parameters are the radius of thevirus, the receptor density ratio and the mobility of the receptors. Figure 15shows the influence which each of these parameters has on the required time.The influence of the radius (Fig. 15a) is analyzed for three different densityratios. The red dashed lines correspond to the maximum radius determinedaccording to Eq. (43). All curves show an increasing trend with the highestvalue for the maximum radius. The curve with a smaller value for ˜ ξ indicatesan increase in time, due to the higher difference in the receptor density availableand required. In Fig. 15b, the time depending on the density ratio for differentradii is shown. A longer process time is required for larger radii. The differencebetween the curves is small for larger ratios and becomes more significant as theratio becomes smaller.Finally, Fig. 15c and 15d both show the influence of the mobility on therequired time. The curves in Fig. 15c corresponds to different density ratios andindicate that lower density ratios are related to the higher time requirements.The difference between the curves for the different densities is more pronouncedas the mobility takes smaller values. A similar behavior is presented in Fig.15d, where the curves correspond to different radii. In both cases a decreasein the mobility causes an increase in the required time, which can be expectedsince a larger m enables the receptors to move more rapidly to the adhesionzone. An analysis of the influence of the radius and of the density ratio to therequired time is also presented in the works by Gao et al. [20]. This research26 a) 0 . . . . . . . T i m e t [ s ] Virus radius R [ µ m]0 .
010 0 .
015 0 .
020 0 .
025 0 .
030 0 . ξ = 1 / . ξ = 1 / . ξ = 1 / . . . . . . T i m e t [ s ] Density ratio ˜ ξ [ − ]0 . . . . . R = 0 . µ m R = 0 . µ m R = 0 . µ m(c) 2 . . . . . . T i m e t [ s ] Mobility m [ µ m / s]0 . . . . . ξ = 1 / . ξ = 1 / . ξ = 1 / . . . . . . . T i m e t [ s ] Mobility m [ µ m / s]0 . . . . . R = 0 . µ m R = 0 . µ m R = 0 . µ m Figure 15: (a) Duration of the process depending on the virus radius for differentreceptor density ratios ˜ ξ . Vertical dashed lines indicate the maximum radii forparticular critical receptor density ratios. (b) Duration of the process dependenton the receptor density ratio ˜ ξ for different virus size R . (c) Process durationversus mobility parameter m for different receptor density ratios ˜ ξ . (d) Processduration versus mobility m for different virus size R .27 a) 0 . . . . . M a x i m u m v i r u s r a d i u s R m a x [ µ m ] Binding range l we [nm]0 . . . . ξ eq req = 7200 [ µ m − ] ξ eq req = 4800 [ µ m − ] ξ eq req = 3600 [ µ m − ] (b) 0 . . . . . T i m e t [ s ] Binding range l we [nm]0 . . . . ξ eq req = 7200 [ µ m − ] ξ eq req = 4800 [ µ m − ] ξ eq req = 3600 [ µ m − ] Figure 16: (a) Maximal virus radius depending on the binding range l we forgiven required receptor density ξ eq req . (b) Process duration depending on thebinding range l we for given required receptor density ξ eq req .group also shows a strong increase in the required time for an increasing radiusbut proposes a minimum value for the radius. Furthermore, contribution [20]shows a rapid increase in the required time when the radius comes close to theminimum value, whereas the current model does not predict such a behavior.Similar to the results which are shown in Fig. 15b, Gao et al. [20] observe astrong decrease in the required time for decreasing density ratios with an upperlimit at ˜ ξ = 1. A more comprehensive study of the uptake process also requires the data oncooperativity (Sec. 8) to be included in the model. The influence of this factoris demonstrated on the basis of two examples dealing with the effects of thebinding range l we . The required receptor density Eq. (42) in combination withthe equation for the radius Eq. (43) provides a relation between the bindingrange and the limiting radius. The corresponding results are presented in Fig.16a where the single curves are related to different receptor densities of the virus.As already shown in Fig. 7b, a smaller receptor density of the virus benefits theprocess, resulting in a larger possible radius.Furthermore, the influence of the binding range l we on the required density(42) and indirectly on the duration of the process is presented in Fig. 16b.Again, the three curves correspond to different receptor densities of the virus.For the upper curve, corresponding to ξ eq = 7200 µ m − , a large binding rangehinders the begin of the process, since the required density ξ eq req becomes to betoo high. On the contrary, the lower curve corresponding to ξ eq = 3600 µ m − ,shows a lower threshold for the process to take place. Here, the duration of the28 a) T i m e t [ s ] Cell density ξ [ µ m − ] Mobility m [ µ m / s] . . . . . . . . . . . . (b) T i m e t [ s ] Mobility m [ µ m / s] Binding range l we [nm] . . . . . . . . . . . . . . . . . Figure 17: (a) Process duration depending on mobility m and the initial receptordensity ξ . (b) Process duration depending on mobility m and the binding range l we . Plots also show the admissibility ranges for chosen parameter sets.process strongly depends on the binding range itself and prefers lower values inorder to complete the process quickly. However, for an extremely small l we , thenumber of virus and cell receptors coming into contact is not sufficient and theprocess cannot start. The discussion of results closes by presenting the interaction of selected processparameters and their influence on the initiation and duration of the process.Figure 17a shows the combination of the initial receptor density ξ and mobility m , and shows a strong increase in the required time for the parameters cho-sen. If both parameters take small values, the process does not start, but oncethis threshold is surpassed the required time drops rapidly regardless of whichparameter is changed. Especially the area with small values for the mobilityis interesting, since it is not uncommon for viruses to attach to cell receptorswith small mobilities. Some typical examples are the HIV-virus connecting to areceptor with mobility m = 0 . µ m / s [47] or the Semliki Forest virus connect-ing to a receptor with m = 0 . µ m / s [20]. Finally, Fig. 17b shows a range ofcombinations for mobility m and binding range l we in which the process takesplace. The influence of the binding range on the required time is weak, com-pared to the influence of the mobility, and becomes more noticeable for smallervalues of the mobility.
10 Conclusion and outlook
The present study focuses on the investigation of the viral entry driven by thereceptor diffusion using the finite difference method as simulation technique. Anapproach based on the consideration of the energetic aspects yields a formulationproviding a well-posed description of the endocytosis process. The motion of the29eceptors is described by the diffusion differential equation accompanied by twoboundary conditions dealing with the flux balance at the ends of the consideredarea. In addition, a supplementary condition is introduced to define the energybalance at the adhesion front.The model developed shows several important features: the definition of thesupplementary condition only depends on the quantities at the front, and thenumerical simulation of the problem bypasses the introduction of assumptionstypical of an analytical solution. The approach is highly efficient with regard totime and computer capacity, such that a fast simulation of different scenarios anda profound study of process parameters are possible. Here, the influences and theinteraction of mobility, receptor densities, virus size and receptor cooperativityplay a central role. Their analysis, for example, yields data on the admissibleregions, the upper limit of the size of the virus able to enter the cell and theestimation of the process duration. Amongst others, the analysis shows thatthe process duration strongly increases when a virus size approaches a criticalvalue and that extremly high and low values of binding range have an impedinginfluence on the process initiation.The results presented in this work pertain to a helical and a spherical viruspenetrating a flat cell surface, which enables the taking of advantage of theaxial and rotational symmetry and perform simulations in a two dimensionalsetup. However, an extension to a three dimensional setup has to be taken intoaccount in order to analyze the receptor distribution for a non-spherical virusor a non-homogeneous receptor density of the cell. Furthermore, additionalcontributions, for example, caused by bending of the cell ahead of the front, canbe considered in the energetic supplementary condition. Alternative expressionsfor bending lipid bilayers can also be introduced in order to carry out morerealistic simulations.
Acknowledgments
We gratefully acknowledge the financial support of the German Research Foun-dation (DFG), research grant No. KL 2678/7-1, and the Austrian Science Fund(FWF), research grant No. I 3431-N32. We also thank Dr. Matias Hernandezfrom the Max Planck Institute of Molecular Physiology at Dortmund, Germany,for valuable discussions.