Near the jamming transition of elastic active cells: A sharp-interface approach
Yony Bresler Benoit Palmieri Martin Grant





=-1











ifundefinedtextcolor

 
 
 
 
 
 
 
 












@breqn@environ#1#2
    breqn@oldbegin@#1#1breqn@oldend@#1end#1breqn@begin@#1
        #1breqn@oldbegin@#1
        \begin#2
    
    breqn@end@#1
        end#1breqn@oldend@#1
        \end#2
        #1breqn@begin@#1
        end#1breqn@end@#1
    
    #1breqn@begin@#1end#1breqn@end@#1
@breqn@environequationdmath
@breqn@environequation*dmath*

    martin.grant@mcgill.caPhysics Department McGill University Montreal Canada
        We use a sharp interface model for active cells to study the jamming
        transition point and behavior near it by varying cell concentration,
        active velocity and elasticity, including a binary mixture of soft
        and stiff cells. We determine the jamming transition point, as well
        as behavior near the transition, including the effective diffusion,
        and sixfold bond correlations. Finally, we expand on previous studies
        by showing the Voronoi dimensionless cell shape can be treated as
        an order parameter at any concentration.
        
          Usage Secondary publications and information retrieval purposes.
          PACS numbers May be entered using the environment .64.70.pm, 87.18.FxNear the jamming transition of elastic active cells: A sharp-interface
        approach
            Martin Grant
            December 30, 2023
        ===============================================================================
        
        
        Assemblies of cells can behave as a liquid, or as a solid. Recently,
        the study of how biological processes relate to concepts such as jamming
        and the glass transition have gained attention. These studies include
        experiments, such as jamming during cancer invasion <cit.>,
        ongoing jamming in the corneal endothelium <cit.>, self-driven
        jamming in microbial growth <cit.>, and unjamming in
        wounds <cit.>. There have also been theoretical biophysics
        studies, which have established a connection between cell migration
        and glass dynamics <cit.>, direct measurement of the
        transition with varying concentration for soft disks <cit.>,
        and a proposed jamming phase diagram <cit.>. Notably,
        the Vertex <cit.> and Self-propelled Voronoi <cit.>
        models showed that cell shape can be used to determine the jamming
        state. These studies have also recently been extended to 3D <cit.>. 
        
        Though driven and active system are far from equilibrium, they have
        been shown to share key features with the equilibrium glass transition,
        such as effective thermal behavior and time correlation functions
        <cit.>. However, despite a wide range of
        experimental and theoretical work on the glass transition, several
        fundamental questions remain unanswered or contested <cit.>.
        In two dimensions, one must consider a potential Kosterlitz-Thouless-Halperin-Nelson-Young
        hexatic phase <cit.>. Indeed, this
        has been seen, for example in colloidal films <cit.> and
        in repulsive Monte Carlo simulations <cit.>, although
        for hard-core or very short-range interactions, transitions appear
        to be first-order <cit.>. Whatever the case may be for
        active systems, these results for thermal systems will play a role.
        
        In this letter we wish to extend these studies using our recently
        developed model for elastic cells <cit.>, a sharp interface
        limit of a phase field model of cells <cit.>, which
        shares some similarities with the recently released CellSim3D
        model <cit.>. Points along a cell interface 𝐑(θ),
        with elasticity γ evolve by minimizing the local curvature
        K, setting a preferred cell area π R_0^2, in the presence
        of repulsion forces due to the other cells, and the cell velocity.
        More precisely, the time evolution for cell n is given in dimensionless
        units by
        
        
            ∂_t𝐑_n(θ,t)=[γ K(θ,t)-1/R_0+μ(A_n-π R_0^2)+β(tanh(α d_n,θ)-1)^2]𝐧_θ+𝐯_n,
        
         where K(θ,t) is the local curvature of the interface,
        1/R_0 is the natural curvature of a spherical cell, μ
        sets the strength of the preferred cell area, A_n is the current
        area, λ is the interface thickness, α=√(15/2)/λ,
        β=150/αλ and d_n,θ is the distance
        to the nearest neighboring cell along the local normal to the interface
        𝐧_θ. The cell velocity 𝐯_n is the
        sum of active and inactive terms: The active velocity has a constant
        magnitude v_A, and reorients with probability P(t)=1/τe^-t/τ,
        where τ is the mean time between reorientations. Inactive velocity
        is due to forces exerted by the other cells surrounding it as well
        as the substrate and surrounding water. 
        
         The model can be used at any degree of confluence, or concentration,
        ρ=Nπ R_0^2/L^2, where N is the number of
        cells and L is the length of the simulation box. This does not
        take into account the distance between adjacent interfaces which is
        proportional to λ, and hence it under-represents slightly
        the actual degree of confluence. Cell stiffness was increased three-fold
        from our previous study, γ=1.35 for 'soft' cells and γ=3.75
        for 'hard' cells, while v_A was varied throughout. The rest of
        the simulation parameters were unchanged, λ=7, R_0=λ^2, μ=0.5, τ=10^4, ϵ=10^3
        and time integration step dt=0.1, all given in a.u.. Each simulation
        consists of 288 cells in a square box with periodic boundary conditions.
        Simulations were run for at least t=800,000 after an equilibration
        of t=80,000 from an initial hexagonal configuration. We expect
        that our results extend beyond the exact details of the model as the
        long-time behavior has been shown to be independent of model details
        in a dense liquid near the glass transition <cit.>.
        
        We begin by demonstrating that our system can reach the jammed state.
        Under the right conditions, the active velocity is unable to allow
        a cell to squeeze through its neighbors leading to dynamical arrest.
        Fig. <ref> shows a snapshot of
        cell interfaces and their center of mass movements over t=8·10^4,
        showing most cells do not exchange neighbors. Of course, that alone
        is not sufficient evidence of jamming. As our system is active rather
        than thermal, we can cross the jamming point by varying v_A,
        the active motor strength. Keeping other parameters constant, Fig.
        <ref> depicts the mean square displacement
        (MSD) over time as a function of simulation time. For short times
        t≪τ, all simulations show ballistic diffusion where MSD  v_A^2t^2.
        By t≫τ however, the behavior depends on v_A. Large active
        velocity (shades of green) have the expected diffusive behavior, MSD ∝ t,
        and the magnitude decreases with v_A. Eventually v_A drops
        below a threshold and MSD becomes sub-linear (shades of red) indicating
        possible dynamical arrest due to caging effects. 
        
        A more direct signature of a transition from a liquid to a jammed
        state is in the falloff of the bond correlations of the sixfold orientation
        order parameter <cit.>Ψ_6(𝐫)=∑_j,ke^6iθ_jk,
        where the sum is over every neighbor k within 2.5R_0 of cell
        j, and θ_jk(𝐫) is the angle between
        them relative to the x-axis. A solid exhibits long range orientational
        order, and the bond-orientational correlation function g_6(r)=⟨ψ_6(𝐫)ψ_6^*(0)⟩
        approaches a nonzero constant for large r, whereas in a liquid
        it decays exponentially with a characteristic correlation length ξ,
        g_6(r)∼ e^-r/ξ. However, in between these may
        lie an additional hexatic phase, where bond orientation decays algebraically,
        g_6(r)∝ r^-η(T) . The Kosterlitz-Thouless
        theory also predicts the hexatic-liquid transition occurs when η=1/4
        (and implicitly, the solid-hexatic transition is at η=0) <cit.>.
        Figure <ref> shows the bond orientational
        correlation function for a series of simulations. Low v_A are
        below the hexatic liquid transition shown as the dashed line. η
        grows with increasing motor strength, crossing into the liquid phase
        where the exponential fall-off would be faster than linear on this
        log-log plot. Given our system size, it is not possible to discern
        whether there exists an additional solid-hexatic phase transition.
        
        We determine the transition point by examining the average effective
        self-diffusion constant, D, as a function of v_A in Fig. <ref>.
        It is computed using the Green-Kubo relation D(t)=1/2∫_0^t⟨ 𝐯(𝐭')· 𝐯(0)⟩𝐝𝐭'
        and averaging over late times. We see D decrease linearly with
        v_A, until it vanishes below a noise floor of D<2·10^-3.
        Using a linear best fit we label the point at which diffusion vanishes
        as v^*. We will show that this coincides with other signatures
        of a jamming transition. See supplemental information for simulations
        videos below (v_A=0.008), very near (v_A=0.009) and above
        (v_A=0.010) the jamming transition. Note these include the equilibration
        time, which is discarded from all statistics.
        
        We use a similar idea to the Turnbull free volume model to obtain
        the effective diffusion for a thermal liquid near jamming <cit.>:
        consider a single cell with velocity v_A moving towards a fixed
        channel of width d<2R+λ. For the cell to squeeze through,
        the dominant term of Eq. <ref> that must be overcome
        is due to the curvature of the cell. Hence, there exists a minimum
        velocity needed to fit through, which is given by v^*∝λ/d.
        If v_A≤ v^*, the cell will not be able to squeeze into the
        channel and become jammed. Higher motor velocities will be able to
        squeeze through and have net velocity v_eff=v_A-v^*, and
        indeed, we see this in simulations of an isolated cells and immovable
        square channel. Given this effective velocity, and that the system
        is still diffusive, we obtain D∝(v-v^*). As
        we will show, this model is also effective at predicting the transition
        velocity as well as diffusion in the dilute regime.
        
        Another method for detecting the transition expands on a result from
        our previous studies <cit.> by investigating the effective
        diffusion for a fixed v_A as concentration is varied. Fig. <ref>
        shows the result remains linear, but now we must account for jamming,
        and 
        
            D(ρ)=Max{ D_0(1-ρ/ρ^*),0} ,
        
         where D_0=1/2v_A^2τ is the diffusion of an isolated
        cell, and ρ^* is the concentration at which jamming occurs,
        i.e. v_A=v^*. The single parameter lines of best-fit show good
        agreement and suggest a method for estimating the onset of jamming
        purely from the diffusion of unjammed cells in a more dilute concentration.
        This differs from previous results of self-propelled soft disks <cit.>,
        we expect due to details in the self-propulsion model.
        
        We performed simulations over a range of concentrations (0.92%-0.80%)
        and mixing ratios of soft to hard cells (0,25,50,75 & 100%).
        For each parameter set, v_A was varied and a linear fit of D
        for values above the noise floor were used to estimate the transition
        point, where D(v^*)=0. The results are summarized
        in Figure <ref>. Each series
        represents different mixing ratios of soft and hard cells. As expected,
        both higher concentration and harder cells lead to transition at higher
        velocity. This reaffirms some of the known behavior of the jamming
        phase diagram <cit.>. Using the same simple model of the
        cell squeezing through a channel, and the free volume between adjacent
        hexagonal cells, we obtain
        
            v^*(ρ)=γ̅2·3^1/4/√(π)(R+λ/2)(√(ρ)-√(ρ_J)),
        
        where γ̅ is the average elasticity of all cells, and
        ρ_J is the concentration of point J where jamming
        occurs at zero velocity. This single parameter fit is shown as a dashed
        line and is in good agreement with our results, particularly for the
        harder cells and mixed cells, and gives similar estimates for ρ_J.
        This is within error of previous results <cit.>
        for soft repulsive discs, ρ_J=0.843, since our definition
        of ρ does not include the interface width. The soft cells, as
        well as the results at the lowest concentration, appear to have larger
        systematic error, likely due to the smaller values of v* approaching
        our noise floor. 
        
        Having determined v^* for each series, we now examine behavior
        across the jamming transition, v_A-v^* in Fig. <ref>.
        In all plots, symbol color denotes concentration, while the symbol
        shape denotes hard to soft cell mixing ratio, as shown in the legend.
        Figure <ref> shows a compilation for the
        effective diffusion constant. D
        increases with a universal slope above the transition, as seen by
        the line of best fit. Once again, this is consistent with the prediction
        of our cell squeezing model. Furthermore, the best-fit slope D(v-v^*)=8.264(v-v^*),
        is surprisingly close to half the Turnbull coefficient <cit.>,
        2R/6≈16.33. In the case where both soft and hard cells are
        present, the average value is shown. The inset shows the difference
        in diffusion constant between soft and normal cells of the same simulation,
        where the y-axis is 10 times smaller than the main figure in order
        to highlight the small differences. Although there might be a slight
        bias for softer cells to have higher diffusion than normal cells in
        the same simulation at large v_A-v^*, the cells appear to unjam
        simultaneously as there is no discernible difference near v_A-v^*≈0. 
        
        Fig. <ref> shows a similar compilation for η.
        Here we find similar but non-universal scaling. Our results appear
        consistent with the KT prediction η(v-v^*=0)=1/4,
        though due to systematic error in determining v^* and η,
        we can not concretely rule out other results such as η=0.385
        which has been suggested in other work <cit.>. As η
        is poorly defined in the liquid state, the inset shows that the above
        v_A-v^* the liquid correlation length ξ combined with
        the average elasticity appears to be universal, as γ̅/ξ^2∝ v_A-v^*. 
        
        Bi et al. showed <cit.> that using a self-propelled
        Voronoi, the average dimensionless shape index, ⟨ q⟩ =⟨ P⟩ /√(⟨ A⟩)
        (where P is the cell perimeter) was a universal order parameter.
        They found that at and below the transition exists a universal value
        ⟨ q⟩ _SPP=3.81 , which increases above
        the transition. We see that the same does not hold for this model.
        The two panels show the average cell shape for hard (Fig. <ref>)
        and soft (Fig. <ref>) cells. Results below the
        jamming transitions are sensitive to equilibration procedure hence
        they are shown with hollow symbols. The figures show that while the
        shape index is increasing above the transition, the value varies significantly
        between soft and hard cells, as well as with concentration and mixing
        ratio, and hence is not a unique order parameter. 
        
        Of course, the self-propelled Voronoi model differs from our sharp
        interface in several ways, most significantly simulations are always
        performed at full confluence ρ=1, cells have limited deformation,
        and all cells had the same elasticity. We make a more direct comparison
        by computing the Voronoi tessellation of the cell center of mass from
        our simulation, similar to that used to analyze the local structure
        and statistics of experimental cell data <cit.>. As shown
        in Fig. <ref>, the Voronoi shape index does appear
        to be a good order parameter for the transition. But our value at
        the transition q_V^*≈3.765 is lower than the self-propelled
        Voronoi model value (q_SPP^*≈3.81, shown as dashed line).
        Note that for v_A-v^*<0 the shape index appears to go below
        the value at the transition for some simulations, in contrast to the
        self-propelled Voronoi model. We suggest that although the Voronoi
        shape index may be a useful order parameter, the transition value
        appears to be model dependent, and hence may vary between cell lines
        as well.
        
        In conclusion, we have demonstrated that our elastic cell model reproduces
        both liquid and hexatic phase. For each parameter set, we estimate
        the velocity at the onset of jamming v^*, and our results are
        consistent with the hallmarks of jamming. Using a simple free volume
        argument along with simulation results, we have also shown the following
        relations: an exact relation for the transition velocity, shown in
        Eq. <ref>. That the effective diffusion is zero in the
        hexatic phase and has universal linear behavior with increasing v_A,
        and that that diffusion for a fixed v_A below the jamming concentration
        depends on the jamming value. The hexatic-liquid transition appears
        consistent with KT theory η=1/4, and in the liquid phase
        γ̅/ξ^2∝ v_A-v^*. Finally, we showed the
        shape index q of the actual cells (q_S,q_H) at
        the transition varies with system parameters. Using the shape index
        of the Voronoi tessellation, the transition appears constant at all
        concentrations and mixing ratios, q_V^*≈3.765, though
        this is lower than the previously predicted universal value.
        
        We gratefully acknowledged The Natural Sciences and Engineering Research
        Council of Canada and the Fonds québécois de la recherche
        sur la nature et les technologies for funding this research, as well
        as Calcul Québec and Compute Canada for providing computing
        facilities.
        
        aipauth4-1
