Elasticity tunes mechanical stress localization around active topological defects

Abstract

Mechanical stresses are increasingly found to be associated with various biological functionalities. At the same time, topological defects are being identified across a diverse range of biological systems and are points of localized mechanical stress. It is therefore important to ask how mechanical stress localization around topological defects is controlled. Here, we use continuum simulations of nonequilibrium, fluctuating and active nematics to explore the patterns of stress localization, as well as their extent and intensity around topological defects. We find that by increasing the orientational elasticity of the material, the isotropic stress pattern around topological defects is changed substantially, from a stress dipole characterized by symmetric compression–tension regions around the core of the defect, to a localized stress monopole at the defect position. Moreover, we show that elastic anisotropy alters the extent and intensity of the stresses, and can result in the dominance of tension or compression around defects. Finally, including both nonequilibrium fluctuations and active stress generation, we find that the elastic constant tunes the relative effect of each, leading to the flipping of tension and compression regions around topological defects. This flipping of the tension–compression regions only by changing the elastic constant presents an interesting, simple, way of switching the dynamic behavior in active matter by changing a passive material property. We expect these findings to motivate further exploration tuning stresses in active biological materials by varying material properties of the constituent units.

---

1 Introduction

Mechanics have been shown to play an indispensable role in governing cell behavior, from differentiation of stem cells,^1,2 to the migration of epithelial cells,³ and the establishment of bacterial size and shape.⁴ Mechanical stimuli and cell response are connected through the process of mechanotransduction, which translates physical forces into biochemical responses in eukaryotes and prokaryotes.^5,6 Such biochemical responses, in turn, activate signaling pathways leading to various feedbacks to the cell behavior, controlling vital cell functions such as migration, proliferation, and cell death.^7,8 Mechanotransduction plays a role at subcellular scales, controlling the behavior of the cell cytoskeleton,⁹ and at the tissue scale impacting collective cell behaviors such as jamming,¹⁰ wound healing¹¹ and cell extrusion.^3,12 As such, identifying regions of mechanical stress localization within living materials is of considerable importance in various biophysical applications, as these localized stress areas present hot spots of signaling activation in a wide range of living materials.

Recently, through the analogy between several distinct living materials and liquid crystals, topological defects have emerged as regions of stress localization with potential biological functions.¹³ Topological defects – singularities in the orientation field – have been found to be at the core of several important biological functions,¹⁴ such as cell death and extrusion,¹³ mound formation in stem cells,¹⁵ bacterial competition,¹⁶ limb origination in morphogenesis,¹⁷ and differentiation of myoblasts,¹⁸ among others.

In two-dimensional nematic systems, topological defects are generated as pairs of ±1/2 charged defects.¹⁹ As regions of high deformation, defects allow for local phenomena, not possible in the strongly aligned nematic phase. In active systems, while the trefoil-shaped −1/2 defect moves diffusively, the +1/2 has a polar, comet shape, which allows the defect to move along its head–tail axis,²⁰ where the ‘head’ consists mostly of bend deformation and the ‘tail’ consists mostly of splay deformation. The isotropic stress of the +1/2 defect has a distinctive dipole pattern,¹³ generating regions of compression and tension around the defect core. Therefore, the defect – a feature of the orientation field of the cells – by way of generating flow and stress features, affects the mechanics of the cell layer in a highly local way and becomes an actor in the mechanotransduction picture of tissue biology.

This is most readily demonstrated in the case of cell extrusion in mammalian epithelia (MDCK cells),¹³ where +1/2 topological defects with a dipolar tension–compression stress pattern were shown to cause extrusion of cells from the layer. Likewise, tension at the center of a −1/2 defect was shown to lead to hole formation in MDCK monolayers.²¹ More generally, in the tissue bulk, compression has been widely linked to cell removal, while tension has been associated with cell division.²² While significant focus has been put on identifying topological defects and their potential roles across biological matter, much less is known about what controls the localized stress patterns around the defects. In particular, to connect the topological defects to biological functionality, it is of considerable importance to understand how the intensity of the stress and its spread around defects are affected by the intrinsic properties of the living material such as its elasticity. A study of the stresses generated around topological defects is therefore in order.

In passive liquid crystals, orientational elasticity – the resistance of the material to orientational deformations – is known to play a central role in the generation of both flow and stress of topological defects.^23–25 Moreover, in most experimental realizations of passive liquid crystals, the elastic response to orientational deformations is shown to be characterized by two distinct elastic constants that penalize bend and splay deformation modes.¹⁹ The difference in splay and bend elastic constants is especially relevant around +1/2 topological defects because the defect combines a bend region at the head and a splay region at the tail. Indeed, experiments on actin filament-myosin motor protein mixtures have shown that increasing the bend elastic constant by adding microtubules results in alteration of the +1/2 defect shape from a rounder U shape for a high splay elastic constant to a pointed Λ shape for a high bend elastic constant.^26,27 A similar change in defect shape has been observed in continuum and particle-based simulations.^27,28 Nevertheless, the majority of theoretical and modeling studies have been limited to using the single elastic constant approximation^20,24,29,30 and important questions about the role of elasticity on governing the mechanical stresses around topological defects remain unanswered.

Here, we focus on addressing these questions using a continuum model of a nematic liquid crystal. In particular, we explore the impact of orientational elasticity on patterns of stress localization around topological defects in the presence of nonequilibrium fluctuations and active stresses. In both cases, we find a significant increase in the intensity and spread of mechanical stress around defects upon increasing elasticity and show an intriguing crossover from a stress dipole characterized by compression–tension regions around the defects, to a localized stress monopole at the defect position. Furthermore, going beyond the single elastic constant approximation, we find that different bend and splay elastic constants can result in a substantial alteration of stress patterns around topological defects leading to a switch from intense tension and diffuse compression to intense compression and diffuse tension. Combining both active fluctuations and active stresses, we find that the elasticity can tune the relative effect of each, making it possible to invert the stress pattern with elasticity alone.

The paper is organized as follows: first we review the continuum nematohydrodynamics equations and the simulation method (Section 2). This is then followed by the results of varying the single elastic constant for a fluctuating nematics in Section 3.1. We then present results of varying bend and splay elasticities in Section 3.2. Finally, we provide characterization of elasticity effects in the presence of active stresses in Section 3.3, show the effect of varying the single elastic constant on −1/2 defects in Section 3.4, and discuss the relevance of our findings and potential experimental tests of our predictions in Section 4.

2 Methods

We implement fluctuating two-dimensional nematohydrodynamics as described in ref. 25. The equations are implemented with a hybrid lattice-Boltzmann approach.³¹

The nematic order parameter is described by a tensor where n_i is the director, q is the magnitude of the nematic order and δ_ij is the Kronecker delta. The nematic order parameter evolves according to the Beris–Edwards equation:³²

where γ is the rotational diffusivity, andis the molecular field, describing the relaxation of the order parameter to the minimum of the free energy. The superscript ST denotes the symmetric, traceless part. The free energy is comprised of a Landau-de Gennes term and an elastic term:

The Landau-de Gennes term follows from the expansion of the order parameter, Q_ij:

f_LdG = A(1 − Q_ijQ_ji)², (4)

where A is the bulk free energy strength. This constant is chosen such that the nematic state is favoured. The elastic free energy can be obtained by expanding in derivatives in Q_ij as follows:^32,33The well-known Frank elastic energy, with splay, twist, and bend elastic constants k₁, k₂, k₃:¹⁹can be recast in terms of L_i from (5) as:^26,32,34,35When k₁ = k₂ = k₃, we recover the commonly used one constant approximation with elastic constant K:^20,36

Finally, on the left-hand side of the order parameter evolution (1), the co-rotation term S_ij describes the interaction between Q_ij and gradients in the flow v_i:

with the rate of strain tensor E_ij = 1/2(∂_iv_j + ∂_jv_i) and vorticity tensor Ω_ij = 1/2(∂_iv_j − ∂_jv_i). The flow alignment parameter λ controls the alignment of the nematic with the flow gradients, and therefore has an important role in determining the type of dynamics that arise.^37–40

The velocity field is determined by the incompressible Navier–Stokes equations, with density ρ, velocity v_i and generalized stress Π_ij:

ρ(∂_t + v_k∂_k)v_i = ∂_jΠ_ij, ∂_kv_k = 0. (10)

The stress, Π_ij, includes a pressure term, Π^p_ij = −pδ_ij, viscous stresses, Π^v_ij = − 2ηE_ij with viscosity, η, and elastic stresses which is a function of the flow alignment parameter, λ and the free energy density f. In an active nematic, the active term Π^a_ij = −ζQ_ij with activity, ζ, is added. When the activity is positive, bend perturbations in the nematic state are unstable and lead to an extensile active turbulence.²⁰ For a negative activity, splay perturbations are unstable and a contractile active turbulence develops. It is therefore expected that the response of an active nematic to varying bend and splay elastic constants will depend on the strength and sign of the active stress. To isolate the effect of elasticity from that of active stress, we first focus on the breaking of detailed balance by nonequilibrium fluctuations, and then in Section 3.3 we study the synergistic effects of elasticity and active stress. Fig. 1 shows representative snapshots of the active turbulence and defect density evolution with time, comparing systems driven out of equilibrium by fluctuations and by active stress.

Fig. 1 Defect turbulence phases for fluctuating and active nematics. (A) and (B) Show the defect turbulence phase for fluctuating and active nematics respectively, with directors in black, +1/2 defects as green circles, −1/2 defects as blue triangles, and the isotropic stress σ_iso as a heat map. (C) and (D) Show defect density for fluctuating and active nematics respectively, as a function of simulation time nondimensionalised by the passive nematic timescale τ_p = γ/A. (A) and (C) Are performed at ξ_θ = 0.01. (B) and (D) Are performed at [small zeta, Greek, tilde] = 12.5 and both at a representative elastic constant K̃ = 2.5. The results in C and D are averages over 10 replicates and error bars correspond to standard deviations.

2.1 Fluctuations

It has been shown that breaking detailed balance by fluctuations of different kinds, such as fluctuations in the order parameter field, velocity field or angle of the director, can lead to very similar defect dynamics as observed with active stresses.²⁵

Therefore, to break detailed balance without including active stresses, we include fluctuations in the angle of the nematic. At every step we calculate the angle θ_l of the director at lattice point l, and add a random angle, taken from the uniform distribution U, scaled by the fluctuation strength ξ_θ:²⁵

and then reconstruct Q_ij. This method preserves the symmetry and tracelessness of the order parameter, changing only the effective angle.

2.2 Simulations setup

We numerically solve the governing equations over a square with side length L with periodic boundary conditions. The system is initialized with a velocity field at rest and a small amount of noise in Q_ij, as in eqn (11) but with θ(t = 0) = 0. The parameters and their values are listed in Table 1. Simulations were performed over T_sim = 3 × 10⁴Δt time steps, with a discarded warm up period of T_warmup = 1 × 10⁴Δt. The parameters are chosen such that the system is at very low Reynolds number (Re ∼ 10⁻⁴) and the nematic coherence length scale is much smaller than the system size (l_Q/L < 10⁻²). Unless otherwise stated, we report simulation results in dimensionless units. To this end, as the main focus of our quantitative analyses is on the extent (area of influence) of the stresses and their intensity, we define characteristic length and characteristic stress as σ_c = ξ_θA/γ, respectively. Note that since our main control parameter in this study is the elastic constant the characteristic length and stress scales are defined independent of K. As such, the main measurable quantities are reported as follows: the area of influence is normalized by l_c², and the stress intensity by σ_c. Similarly the control parameters are reported as follows using the same characteristic length l_c and characteristic stress σ_c scales: the dimensionless elastic constant is defined as K̃ = K/(σ_cl_c²), and the dimensionless activity as [small zeta, Greek, tilde] = ζ/σ_c. In all simulations we keep the bulk free energy coefficient A, constant, therefore changing elasticity K affects the coherence length . However, this coherence length determines the size of the defect core, i.e. the isotropic region at the core of the defect, while the extent of the stress is predominantly dependent on the average distance between the defects, which is determined by the competition of the activity (from nonequilibrium fluctuations or from active stress) and the elasticity. Furthermore, it is well-established that active defects form even if the energy constant A is zero, through the activity-induced order,⁴¹ which lends further support to the energy constant having only a minor effect on the scale of the patterns around the defect and beyond the defect core. We have further verified that doubling the bulk free energy coefficient A has no significant effect on the extent and intensity of isotropic stresses around topological defects.

Table 1 Simulation parameters with name, symbol and value (or range) and dimension where length: L, mass: M and time: T

Parameter	Symbol	Value	Dimension (2D)
Flow-alignment	λ	1	1
Rotational viscosity	γ	20	M/T
Solvent viscosity	η	40/6	M/T
Density	ρ	40	M/L^2
Bulk free energy strength	A	1	M/T^2
Frank elastic constant	K	[0.01,0.5]	ML ^2/T^2
LB relaxation time	τ LB	1	T
Activity	ζ	[−0.05,0.05]	M/T^2
Initial noise in alignment	n 0	0.05	1
Director angle fluctuation	ξ θ	0.04	M/T
Numerical integration time	Δt	1	T
Square domain length	L	256	L

---

3 Results

Our focus is on the patterns of mechanical stress: the intensity and spread of stress around topological defects. To this end, we conduct simulations of active nematics, first in the presence of nonequilibrium fluctuations, and then including active stresses. In both cases, we perform our analyses in the regime where the system is in a statistical steady state, the topological defect density is saturated and the isotropic stresses are similar (Fig. 1). In order to characterize the stress patterns we calculate the isotropic stress, i.e., the trace of the stress tensor : negative isotropic stress characterizes regions under compression and positive isotropic stress demarcates regions that are under tension. The stress tensor Π_ij accounts for contributions from all components of the stress, including pressure, which ensures incompressibility condition is satisfied, viscous stresses, elastic stresses, and active stresses. The choice of using isotropic stress is motivated by biological implications where tension and compression lead to activation–deactivation of mechanotransductive pathways^8,42 and as a result living materials respond differently to compression and tension.²²

3.1 Elasticity effects on stress localization: single elastic constant

We begin by investigating the effect of varying the elastic constant K, staying within the single elastic constant approximation, on the stress patterns in fluctuating nematics. Since elasticity penalizes deformations in the nematic director field, it is expected that the size of the ordered region in the director field around the defect increases with the increasing elasticity.⁴³ The stress generated around the defect, however, is what influences mechanics, and governs the possible role of topological defects in inducing compression–tension in the material.

At low elasticity, a highly local, dipolar stress pattern is observed with a compressive (σ_iso < 0) and a tensile (σ_iso > 0) region at the head and the tail of the comet, respectively. As such, the defect provides a local region for compression and tension, which are known to differently affect activation–deactivation of mechanosensitive pathways,^8,22 with compression being linked to cell death, extrusion,¹³ and differentiation,¹⁸ while tension has been associated with cell division.²² The compression and tension regions are symmetric around the defect core: they have the same size and have similar stress intensities (Fig. 2). As the elasticity is increased, the compression and tension regions respond differently and the symmetry of the isotropic stress between the head and the tail regions is lost: the tension magnitude increases faster than compression. The pattern then changes from a dipolar isotropic stress pattern, with equal peaks for compression and tension at the head and tail of the defect, to a pattern of a large tensile peak, but a wider and weak compressive peak the defect head. Both tensile and compressive regions expand, but while tension localizes towards the defect core, the compression spreads around the head region of the defect.

Fig. 2 Isotropic stress patterns around defects for varying elasticity K, within one-constant approximation. The upper row shows the increase in stress magnitude, growing by 2 orders of magnitude as K increases from K̃ = 2 to K̃ = 25. The lower row shows that, while for lower elasticity the tensile and compressive regions are of similar strength and size, at higher elasticity tension dominates, changing the character of the stress pattern qualitatively. The average is taken over >1000 defects for low, intermediate, and high elasticity values and 5 replicates were made.

To quantify the effect of elasticity on the intensity and the extent of the stresses around defects, we measure the area of influence and average magnitude of the stress, henceforth called intensity, for both the tensile and compressive regions. To calculate the area of influence, we first find the maximum tension and compression values and measure the size of the area over which the magnitude of the isotropic stress drops to 50% of the maximum value. This gives us a measure of how the spread of both compression and tension regions change with varying the elastic constant. Our quantitative measures clearly show that the areas of influence increase for both compression and tension as the elasticity is increased (Fig. 3). At a low value of the elastic constant, both the tension and compression regions have the same area of influence around the defect. Upon increasing the elastic constant, the size of the tensile area steadily increases, while the compression area expands to a much larger extent, to almost five-fold larger than the tensile area (Fig. 3a). However, measuring the intensity of the stress (Fig. 3b) reveals that at higher elasticities the less spread, tensile area has an almost 10-fold larger stress intensity compared to the more spread, compressive region. Together these results show that, even within the single elastic constant approximation, changing elasticity has a significant effect on the patterns of mechanical stress localization around topological defects: while at low elasticities, both intensity and the area of influence of tensile and compressive regions are symmetric around the defect core, increasing the elastic constant breaks this symmetric localization of the stress leading to a less spread out, highly localized, tension at the defect core with a strong intensity, and a more spread out, but less intense compression at the head of the defect.

Fig. 3 Intensity and extent of stresses around defects increase as elasticity is increased. (A) Area of influence of tension (blue) and compression (red) regions around defect. At low elasticity, both peaks are of similar size, but after a threshold K the compressive peak grows faster than the tensile peak, after which growth stabilizes and both peaks have a fixed difference in size. (B) Intensity of stress in the peak region for tension (blue) and compression (red). At lower elasticity, the peaks have the same strength, however as the elasticity is increased, the tensile peak dominates. Stars indicate values of K̃ used in Fig. 2 Simulation results are averaged over 5 replicates. Error bars indicate standard deviation.

3.2 Elasticity effects on stress localization: distinct splay k₁ and bend k₃ elastic constants

So far, we have only explored elasticity effects on stress patterns around defects, within the single elastic constant approximation, where the nematogens have the same resistance to bend and splay deformation modes. As mentioned earlier, this is typically not the case in passive nematics,⁴⁴ and it has also been shown that cytoskeletal filaments and motor protein mixtures have different bend and splay elastic constants.²⁶ Furthermore, because half-integer defects combine bend and splay deformation modes, we expect the distinction between bend and splay elastic constants to be specifically important there. To this end, we next focus on the impact of elastic anisotropy on the patterns of isotropic stress, their intensity, and spread around topological defects.

The bend-to-splay ratio is defined as κ = k₃/k₁. For low values of κ, when the splay elastic constant is higher, the extent of the stress pattern is small, and the tension has higher intensity than the compression (Fig. 4 top row). As the bend-to-splay ratio is increased to intermediate values κ ∼ 1, the tension has still has intensity, reflecting the results from Section 3.1. As the bend-to-splay ratio is increased to κ ≫ 1, we find that the compression is enhanced until it dominates over tension.

Fig. 4 Isotropic stress when changing bend to splay ratio κ. Top row: The normalized isotropic stress fields show a more intense tensile peak at the tail and a weaker compressive peak at the head for low κ. At high κ, this effect inverts and the compression is more intense while the tension is weaker at the tail. (A) At low κ the compressive and tensile peaks are similar, at intermediate κ compression is larger, and as the ratio approaches κ ≈ 20 tension becomes larger. The stars indicate the values of κ at which the flow fields were measured. (B) When splay dominates, the tensile peak is around ≈2 times stronger than the compressive, with that ratio rising as κ increases. However, as the ratio passes κ ≈ 10, the intensity of compression becomes more dominant. Simulation results represent average fields over >1000 defects and 5 replicates were made. Error bars indicate standard deviation.

Similar to Section 3.1, the effect of changing the elastic anisotropy is best understood by quantifying the intensity and the spread of tension and compression around defects as the bend-to-splay elasticity ratio is varied (Fig. 4).

Three regimes are evident: low, intermediate, and high bend-to-splay ratio κ. The intermediate regime, as expected, is close to the single constant approximation, where the tension has a smaller peak with higher intensity, while the compression has a more diffuse lower peak, as described earlier. However, the low and high κ regimes differ substantially. At low κ, the tensile peak is more intense but has a smaller extent, whereas at high κ it is the compressive peak that is more intense but has a smaller area of influence (Fig. 4). This is an inversion of the effects of tension and compression, although their positioning around the defect core stays the same.

It is noteworthy that the largest stress patterns by area appear in the intermediate bend-to-splay ratio, which also has the largest intensity difference, ≈10-fold, whereas at low and high κ, the areas are small and the intensity differences are relatively low, ≈2-fold (Fig. 4).

3.3 Interplay between activity and elasticity

In all the results presented so far, we focused on fluctuation-induced defects and purposefully neglected dipolar active stresses, which are commonly employed in models of active nematics.²⁰ The purpose was to isolate the effect of elasticity from that of the active stress in inducing coherent stress patterns around topological defects. This is especially relevant since previous theoretical and numerical works within the over-damped compressible limit, and in discrete simulations, have shown that active stresses can effectively normalize the elastic constants.²⁸ Having understood the impacts of elasticity, both at the single elastic constant limit and distinct bend-splay modes, we now turn to characterizing the interplay between activity and elasticity in establishing the stress patterns, their region of influence, and intensity around topological defects. Motivated by experiments on eukaryotic cells,^45,46 and actomyosin cytoskeletal filaments,⁴⁷ here we consider the contractile active stress, i.e. a force dipole along the nematogen's direction that is pointing inwards, and therefore contracts along the director.

In Sections 3.1 and 3.2, we took a systematic approach of first varying only the elastic constant within the single elastic constant approximation and then analyzing the impact of elastic anisotropy. The results showed that the main trends associated with the stress intensity and extent are common within the two approaches. Thus, when we consider the interplay with activity we only focus on the single elastic constant approximation, since it is already expected that activity (depending on its sign) can favor one mode (splay or bend) over the other. Fig. 5 shows the isotropic stress around a defect in a contractile active nematic, ζ < 0, with nonequilibrium fluctuations, ξ_θ = 0.04, within the single elastic constant approximation. At low elasticity, the activity effect dominates and, as expected for a contractile active nematic, the compression → tension dipole is aligned along the tail → head direction around the comet-like +1/2 defect. However, as the elasticity is increased, the compressive and tensile regions flip, resembling the stress pattern expected for an extensile active nematic, although the (contractile) active stress was kept unchanged.

Fig. 5 Isotropic stress patterns around defects for varying elasticity K, and a contractile active stress. (A) Shows the increase in stress magnitude. (B) Shows that while for lower elasticity the compression → tension dipole is pointed in the tail → head direction, at higher elasticity the direction of the dipole flips to tension → compression along the tail → head direction. (C) The sign of the stress gradient at the core of the defect changes with increasing stiffness, indicating that tension and compression have switched places. Orange stars indicate points at which the stress fields are shown. [small zeta, Greek, tilde] = −25, and the averages are made over >500 defects and 5 replicates.

A simple explanation is that while the in an active nematic the stress dipole is dominated by activity,⁴⁸ the stress localization for the fluctuation induced defect is controlled by elastic effects.²⁵ As the elastic constant, K, increases, the strength of elastic effects is increased until they overcome the active effects, and the stress dipole flips direction from compression → tension dipole to tension → compression along the tail → head direction around the comet-like +1/2 defect. This switch in the tension–compression regions by only changing the elastic constant presents an interesting, simple, way for switching the dynamic behavior in active matter by changing a passive material property. It is also noteworthy that this flipping of the stress also corresponds to the switch between extensile and contractile behavior of topological defects in active nematics.²⁰ Experimental evidence of such a switch is presented in epithelial layers upon weakening of cell–cell adhesion,⁴⁹ and while various active mechanisms have been proposed for the flipping of extensile-contractile defect behavior,^25,50 these results show that such a flipping can simply occur upon changing of the oreintational elasticity of the nematogens, a passive material property.

To characterize the flipping of tension–compression regions clearly, we measure the stress gradient along the tail → head direction around the defect core where the average is taken perpendicular to the defect axis and the maximum gradient is picked along the defect axis, and define

such that μ < 0, characterizes tension localized at the tail and compression at the head of the defect, while μ > 0 corresponds to the reverse localization. We see the value of μ switch to negative at K̃ ≈ 2.2 (Fig. 5C) although this value is dependent on fluctuation and active stress strength. Having shown the effect of changing elasticity in a nematic with active fluctuations and stresses, we then construct a phase diagram of the system behavior to summarize the effect of elasticity and active stress.

We find that the phase boundary between tension at the head and tension at the tail is determined both by the elasticity and the activity as ∼K/ζ (Fig. 6). The phase can therefore be changed by varying either of these parameters, both of which are experimentally accessible.

Fig. 6 Phase diagram of the stress gradient around defects in the orientational elasticity-activity phase space. Varying either elasticity or activity can lead to flipping of the direction from tension → compression (blue square) dipole to compression → tension along the tail → head direction (red circle) around the comet-like +1/2 defect. Each point is averaged over 5 realizations.

3.4 Elasticity effects on compression–tension patterns around −1/2 defects

Finally, we explore the effect of varying the single elastic constant, K, on the trefoil-like, −1/2 defects. Since the single and multiple elastic constants approaches give qualitatively similar results, here we focus on the single elastic constant case. Fig. 7A shows the isotropic stress around a −1/2 defect in a passive system (ζ = 0) with nonequilibrium fluctuations. When the elastic constant is low, the bend regions are under tension, whereas splay regions experience compression. As the elastic constant increases the area of tension in bend regions also increases. For the contractile system (ζ < 0) with nonequilibrium fluctuations, in turn, the bend regions are tensile (Fig. 7B). This changes when the elastic constant is increased, where the bend regions are now compressive. This stress inversion is reminiscent of our earlier finding of the inversion of the isotropic stress dipole pattern around +1/2 defect upon increasing the elastic constant for a contractile and fluctuating system.

Fig. 7 Isotropic stress patterns around −1/2 defects for varying single elastic constant, K. (A) In the presence of nonequilibrium fluctuations, for lower values of elasticity, the bend regions around −1/2 are under tension. As the elasticity is increased, the area of tension in bend regions increases. (B) In the active contractile system, at low K the bend regions are compressive, whereas at higher K the bend regions are tensile. The averages are made over 5 simulations.

4 Discussion

Topological defects in the alignment of living entities have been found to play a functional role in a range of biological behaviors by providing local environments for enhanced mechanical stresses.^13,14,51 In this work, we have taken a systematic approach of exploring intensity and area of influence of mechanical stresses around topological defects. To disentangle the effect of self-generated active stresses from passive material properties we separately investigated stress patterns in out-of-equilibrium fluctuating nematics and in active nematics. Our results show that in both cases, the orientational elasticity of nematogens plays a significant role in determining the stress localization patterns around defects. Due to the generic presence of fluctuations in active nematic systems, from subcellular fluctuations in focal adhesion⁵² or stress fibers⁵³ to cell level shape fluctuations in MDCK layers⁵⁴ or Drosophila,⁵⁵ we expect that understanding the interplay between activity and fluctuations will play an important role in understanding cell behavior in the future.

First, in a fluctuating system, we showed that increasing elasticity breaks the symmetry between compressive and tensile regions at the head and the tail of the defect, respectively. While the area of influence of the stress and the stress intensity increase for both compression and tension regions, we find that such enhancement occurs disproportionately: upon increasing elasticity, the region of tension localizes close to the defect core and the tension intensity becomes about an order of magnitude stronger than the compression. At the same time, the compressive region spreads over a relatively larger region around the defect compared to tension. In MDCK cells, knockdown of the mechanosensor α-catenin, and therefore weakening cell–cell contacts, was shown to reduce the overall extent of the isotropic stress patterns around defects, which was interpreted as lowering the elastic constant.¹³

Upon relaxing the one-constant approximation and allowing for different bend and splay elasticities, we find that for low and intermediate bend to splay ratios, the stress patterns around defects are characterized by strong, highly local tension and weak, more spread out compression. However, this behavior switches for high bend to splay ratios where local strong compression and wide-spread weak tension regions emerge.

Finally, combining the effects of activity and elasticity, we present the possibility of the modulation of stress intensity, area of the influence of stress, and also flipping of the compression–tension regions. This inversion of compressive and tensile regions was observed in both +1/2 and −1/2 defects. Taken together, the results show that in a nonequilibrium fluctuating nematic system, we saw that changing the elasticity can change the dipole stress pattern to increase the intensity of tension and the area of influence of compression. In an active and nonequilibrium fluctuating nematic system, we found that changing the elasticity can even lead to the flipping of the localization of compression–tension regions around the defects.

The results presented herein can be useful in understanding how stress localization is controlled by the self-organization of the active entities. Understanding the strength and extent of the localization of stresses that collections of elongated active entities are subject to, and how they depend on material properties of the entities themselves, is essential in studying the mechanically-induced biochemical changes in living biological matter. For example, the compressive region of the isotropic stress pattern of a topological defect has been found to lead to cell death and extrusion in mammalian epithelia¹³ and tension in a tissue is associated with cell division²² or gap opening on soft substrates.²¹ Indeed, the existing experiments that have measured the compression–tension patterns around the defects, have shown that molecular perturbation of cell–cell contacts, which inevitably alters the mechanical properties of the tissues, can affect the intensity and extent of isotropic stress localization around defects,¹³ and could even reverse the direction of the motion of defects.⁴⁹

Experimentally, varying the bend to splay constant is possible, as simply varying the length of a flexible nematic such as actin has been shown to change the bend to splay ratio.^26,28 Furthermore, adding carbon nanotubes to a molecular nematic was shown to increase the bend to splay ratio by up to 6 fold,⁵⁶ while in a similar fashion, microtubules were added to an active nematic consisting of actin and myosin to increase the bend constant >2-fold.²⁶ We predict that the change from intense tension to intense compression when varying the bend to splay ratio should affect cellular systems, as tension and compression are linked to differing biological behaviors.²²

We expect our results to help design new experiments aiming at controlling stress localization, and potentially mechanotransduction, by adjusting the elastic properties of the cell layers.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

A. D. acknowledges funding from the Novo Nordisk Foundation (grant No. NNF18SA0035142 and NERD grant No. NNF21OC0068687), Villum Fonden (Grant no. 29476), and the European Union (ERC, PhysCoMeT, 101041418). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. A. A. acknowledges support from the European Unions Horizon Europe research and innovation program under the Marie Sklodowska-Curie grant agreement No. 101063870 (TopCellComm).

Notes and references

1. M.-K. Hayward, J. M. Muncie and V. M. Weaver, Dev. Cell, 2021, 56, 1833–1847.
2. K. H. Vining and D. J. Mooney, Nat. Rev. Mol. Cell Biol., 2017, 18, 728–742.
3. B. Ladoux and R.-M. Mège, Nat. Rev. Mol. Cell Biol., 2017, 18, 743–757.
4. A. Persat, C. D. Nadell, M. K. Kim, F. Ingremeau, A. Siryaporn, K. Drescher, N. S. Wingreen, B. L. Bassler, Z. Gitai and H. A. Stone, Cell, 2015, 161, 988–997.
5. N. Wang, J. Phys. D: Appl. Phys., 2017, 50, 233002.
6. Y. F. Dufrêne and A. Persat, Nat. Rev. Microbiol., 2020, 18, 227–240.
7. A. Pocaterra, P. Romani and S. Dupont, J. Cell Sci., 2020, 133(2) DOI:10.1242/jcs230425.
8. T. Panciera, L. Azzolin, M. Cordenonsi and S. Piccolo, Nat. Rev. Mol. Cell Biol., 2017, 18, 758–770.
9. S. Mathieu and J.-B. Manneville, Curr. Opin. Cell Biol., 2019, 56, 34–44.
10. D. Boocock, T. Hirashima and E. Hannezo, bioRxiv, 2023 DOI:10.1101/2023.03.24.534111.
11. A. Brugués, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Muñoz, G. W. Brodland, B. Ladoux and X. Trepat, Nat. Phys., 2014, 10, 683–690.
12. S. Monfared, G. Ravichandran, J. Andrade and A. Doostmohammadi, eLife, 2023, 12, e82435.
13. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
14. A. Doostmohammadi and B. Ladoux, Trends Cell Biol., 2021, 140–150.
15. K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327–331.
16. O. J. Meacock, A. Doostmohammadi, K. R. Foster, J. M. Yeomans and W. M. Durham, Nat. Phys., 2021, 17, 205–210.
17. Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun and K. Keren, Nat. Phys., 2021, 17, 251–259.
18. P. Guillamat, C. Blanch-Mercader, G. Pernollet, K. Kruse and A. Roux, Nat. Mater., 2022, 21, 588–597.
19. P. G. D. Gennes and J. Prost, The physics of liquid crystals, Clarendon Press; Oxford University Press, Oxford, New York, 2nd edn, 1993.
20. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246.
21. S. Sonam, L. Balasubramaniam, S.-Z. Lin, Y. M. Y. Ivan, I. Pi-Jaumà, C. Jebane, M. Karnat, Y. Toyama, P. Marcq and J. Prost, et al. , Nat. Phys., 2023, 19, 132–141.
22. T. Zulueta-Coarasa and J. Rosenblatt, Curr. Opin. Genet. Dev., 2022, 72, 1–7.
23. G. Tóth, C. Denniston and J. M. Yeomans, Phys. Rev. Lett., 2002, 88, 105504.
24. D. Khoromskaia and G. P. Alexander, New J. Phys., 2017, 19, 103043.
25. L. Bonn, A. Ardaševa, R. Mueller, T. N. Shendruk and A. Doostmohammadi, Phys. Rev. E, 2022, 106, 044706.
26. R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E124–E133.
27. N. Kumar, R. Zhang, J. J. de Pablo and M. L. Gardel, Sci. Adv., 2018, 4, eaat7779.
28. A. Joshi, E. Putzig, A. Baskaran and M. F. Hagan, Soft Matter, 2019, 15, 94–101.
29. S. P. Thampi, R. Golestanian and J. M. Yeomans, EPL, 2014, 105, 18001.
30. D. Khoromskaia and G. Salbreux, eLife, 2023, 12, e75878.
31. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921.
32. A. N. Beris and B. J. Edwards, Thermodynamics of flowing systems: with internal microstructure, Oxford University Press, New York, 1994.
33. B. J. Edwards and A. N. Beris, J. Rheol., 1989, 33, 1189–1193.
34. K. Schiele and S. Trimper, Phys. Status Solidi B, 1983, 118, 267–274.
35. N. Kumar, R. Zhang, S. A. Redford, J. J. D. Pablo and M. L. Gardel, Soft Matter, 2022, 18, 5271–5281.
36. S. Thampi and J. Yeomans, Eur. Phys. J.: Spec. Top., 2016, 225, 651–662.
37. S. A. Edwards and J. M. Yeomans, EPL, 2009, 85, 18008.
38. K. Thijssen, L. Metselaar, J. M. Yeomans and A. Doostmohammadi, Soft Matter, 2020, 16, 2065–2074.
39. S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, Soft Matter, 2019, 15, 1597–1604.
40. S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, Phys. Rev. Lett., 2020, 125, 148002.
41. S. P. Thampi, A. Doostmohammadi, R. Golestanian and J. M. Yeomans, EPL, 2015, 112, 28004.
42. X. Cai, K.-C. Wang and Z. Meng, Front. Cell Dev. Biol., 2021, 9, 673599.
43. M. Kleman and O. D. Lavrentovich, Philos. Mag., 2006, 86, 4117–4137.
44. D. Demus, J. W. Goodby, G. W. Gray, H. W. Spiess and V. Vill, Handbook of liquid crystals,low molecular weight liquid crystals I: calamitic liquid crystals, John Wiley & Sons, 2011, vol. 2A.
45. O. Du Roure, A. Saez, A. Buguin, R. H. Austin, P. Chavrier, P. Siberzan and B. Ladoux, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 2390–2395.
46. G. Duclos, C. Erlenkämper, J.-F. Joanny and P. Silberzan, Nat. Phys., 2017, 13, 58–62.
47. K. Clark, M. Langeslag, C. G. Figdor and F. N. van Leeuwen, Trends Cell Biol., 2007, 17, 178–186.
48. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365.
49. L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. H. N. S. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, A. S. Yap, Y. Toyama, R.-M. Mège, J. M. Yeomans and B. Ladoux, Nat. Mater., 2021, 1–11.
50. A. Killeen, T. Bertrand and C. F. Lee, Phys. Rev. Lett., 2022, 128, 078001.
51. S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti and V. Vitelli, Nat. Rev. Phys., 2022, 4, 380–398.
52. S. Plotnikov, A. Pasapera, B. Sabass and C. Waterman, Cell, 2012, 151, 1513–1527.
53. L. Guolla, M. Bertrand, K. Haase and A. E. Pelling, J. Cell Sci., 2012, 125, 603–613.
54. S. Zehnder, M. Suaris, M. Bellaire and T. Angelini, Biophys. J., 2015, 108, 247–250.
55. M. Olenik, J. Turley, S. Cross, H. Weavers, P. Martin, I. V. Chenchiah and T. B. Liverpool, Phys. Rev. E, 2023, 107, 014403.
56. S. Turlapati, R. K. Khan, P. Ramesh, J. Shamanna, S. Ghosh and N. V. Rao, Liq. Cryst., 2017, 44, 784–797.

---