Progress toward establishing the mechanical properties of animal tissue is limited by both the availability of measurements and the development of quantitative models that can realistically account for both the material properties and the geometry of cells. To tackle the first challenge, we developed a novel technique which is inspired by previous measurements on tissue culture cells where a bendable cantilever of glass was used as a force probe (Desprat et al., 2005). In these experiments, the deflection of the cantilever is used to quantify the force applied by the tip; in combination with the resulting tissue deformation, a force-displacement relationship for the tissue can be established.

To probe the mechanics of individual embryonic cell edges, we needed to create a novel tool that (1) can be inserted into a single embryonic cell, (2) is sufficiently compliant to appreciably bend under the relevant nanonewton-range forces, and (3) is capable of being visualized within the rather opaque embryonic interior. To meet these needs, we developed a technique for manufacturing a probe consisting of a long stiff glass pipette with a soft, fluorescent polymer tip that acts as a short flexible cantilever (Figure 1a). Briefly, a standard injection glass pipette is used as a mold for making a PDMS (polydimethylsiloxane) cantilever. The tip of a glass pipette is filled with an uncured mixture of polymer, crosslinker, and a green variant of BODIPY (boron-dipyrromethene) dye, which enters by capillary action through an opening at the tip. After the polymer is cured, the tip of the glass scaffold is etched away in a solution of hydrofluoric acid to expose the soft PDMS core. This technique is described in further detail in the Materials and methods section and in Figure 1—figure supplement 1.

Figure 1 with 1 supplement see all

Download asset Open asset

With this probe in hand, we began to study the mechanical properties of cellularization stage Drosophila embryos. From fertilization through the earliest several cycles of mitosis, the Drosophila embryo is a syncytium, that is a very large cell with multiple nuclei sharing common cytoplasm. Subsequently, in the process of cellularization, membranes extend from the surface of the embryo and protrude between the nuclei, such that the cytoplasm is for the first time partitioned into individual cells (Campos-Ortega and Hartenstein, 1985). During cellularization, every newly forming cell is in effect open on its basal side, such that the interior of that cell is connected to the yolk compartment residing in the center of the embryo (see Figure 1b for a schematic illustration). Drosophila embryos at this stage have become a key model for studying epithelial mechanical properties (Bambardekar et al., 2015; D’Angelo et al., 2019; Doubrovinski et al., 2017), as it forms a basis for understanding the molecular and physical changes that must occur to drive the large-scale tissue morphogenesis that begins shortly afterward (Blankenship et al., 2006; Kiehart et al., 2017; Martin et al., 2009).

To introduce our probe into a cellularization-stage embryo, we use a pulled glass pipette to cut a slit in the embryo on the opposite surface from which we plan to take measurements (Figure 1b). Importantly, we found that when a fly embryo is cut along one side, the side that is not cut continues to develop completely normally showing no noticeable defects. We insert the probe through this cut, then through the basal hole of a cell on the far side of the embryo, so that the flexible tip of the cantilever is positioned within a single cell. For consistency, we always make the initial incision on the ventral side of the embryo and probe material properties of cells of the dorsal side. Once the probe is inserted, we use a piezo actuator mounted on the microscope stage to control cantilever movement. As the cantilever is translated parallel to the cellular layer, it pushes against the adjacent cell edge, thus exerting a pulling force on that cell (Videos 1 and 2, Figure 1b).

Video 1

Download asset

Video 2

Download asset

Each pulling experiment consists of three phases. First, during the loading phase, the base of the cantilever is translated at a constant velocity of 0.5 µm/s for approximately 60 s. Next, during the holding phase, the cantilever is held in place for approximately 60 s. Finally, in the unloading (or relaxation) phase, the cantilever is rapidly retracted from the interior of the cell, and the cellular layer is then allowed to freely relax. Through this experiment, the bending of the cantilever and the deformation of the cellular layer are imaged for subsequent quantification. For time frames from a typical pulling experiment, see Figure 1c.

To quantitatively examine displacement information from this type of experiment, we plotted the data in two distinct ways. First, we represented the imaging data from each experiment as a kymograph, with the x-axis indicating time, and the y-axis displaying the position of the cell membranes (as measured from the fluorescent signal intersecting a curve drawn through the middle of the cellular layer). The kymograph for a representative pulling experiment is shown in Figure 2a.

Figure 2 with 1 supplement see all

Download asset Open asset

Additionally, we represent the imaging data from each experiment as a deformation profile. For a given time point, we generate a curve showing the displacement of each membrane from its initial position as a function of its initial position (Figure 2b). Plotting several such curves, each representing a different time point of the experiment, together on one graph allows easy visualization of how tissue displacement spreads over time as the tissue is pulled. The deformation profiles for the loading and holding phases of a representative pulling experiment are shown in Figure 2c and c'', respectively. For some analyses, it is also useful to generate normalized deformation profiles, in which the displacement (or y) axis is linearly scaled so that the maximal displacement for each time point equals 1 (Figure 2c'). Importantly, measurements performed on different specimens are highly consistent (Figure 2—figure supplement 1).

One key advantage of the cantilever technique is that applied force can be calculated from cantilever deflection, which itself can be measured directly from the imaging data. This can be done because cantilever deflection is linearly proportional to the force being applied to it for the range of deflections observed. By Newton’s third law, this force is equal in magnitude to the force the cantilever is exerting on the tissue.

To calculate the force, it is necessary to calibrate the probe. The force-deflection relation of each cantilever can be readily calculated using two pieces of information: the geometry of the individual cantilever, and the Young’s (elastic) modulus of the PDMS material from which it is made. To calculate the Young’s modulus of PDMS, we used three complementary techniques: (1) fabricating a macroscopic rod of PDMS to measure its deformation under its own weight (Figure 3—figure supplement 1a-b), (2) pulling the cantilever through a fluid of known viscosity to measure the resulting deflection (Figure 3a–b), and (3) probing the cantilever with an AFM probe (Figure 3—figure supplement 1c-d). The three methods estimate the Young’s modulus of the PDMS comprising the cantilever to be 1.6 MPa, 0.9±0.17 MPa, and 2.1±0.23 MPa, respectively (the first value is based on a single measurement, see Figure 3—figure supplement 1a-b). The details of these measurements, along with the procedure we used to combine this information with individual cantilever geometry to calibrate each probe, are described in the Materials and methods section. Using this information, we calculated the force applied over the time course for five representative experiments (Figure 3c). The average force measured at maximal tissue displacement (at the end of the loading phase) was 11.1±1.2 nN, and fell to 5.0±0.45 nN by the end of the holding phase.

Figure 3 with 1 supplement see all

Download asset Open asset

In order to correctly infer the various force contributions in a mechanical system, it is important to know whether the dynamics are adiabatic or not. An adiabatic process is defined as one in which no heat transfer takes place. Claiming that a process is adiabatic is equivalent to stating that friction is negligible, or alternatively, that the system is at equilibrium at every time point.

During the holding phase, most of the cell edges remain in a constant position, although the green cantilever trace and a few nearby cell edges move slightly upward (Figure 2a and c''). Since most motion ceases as soon as the cantilever stops moving, the tissue must have been largely in mechanical equilibrium at the time the cantilever stopped moving; in other words, most of the tissue was deforming adiabatically during the loading phase.

Further support for this conclusion comes from examining the loading phase more quantitatively. In previous work, we analytically and computationally explored the behavior of thin material sheets with a variety of mechanical properties (Doubrovinski et al., 2017). When these sheets were pulled at a single point, similarly generated deformation profiles from different time points of the same experiment were invariant up to re-scaling the y- and the x-axis by the factors $t^{\alpha }$ and $t^{\beta }$ respectively (see Figure S4 of Doubrovinski et al., 2017), where the values $\alpha$ and $\beta$ are determined by the mechanical properties of the system. In particular $\beta$ should exceed 0.5 if friction contributes significantly to the dynamics. In the pulling experiments presented here, we find that the deformation profiles obtained from the loading phase (Figure 2c) can be collapsed on almost the same master curve by re-scaling or normalizing the y-axis only; that is, $t^{\beta }=1$, or equivalently $\beta =0$ (Figure 2c'). This indicates that the tissue is deforming without significant input from external friction, which is consistent with the tissue deforming adiabatically (Figure 2d).

So far, we have demonstrated that the tissue is deforming largely adiabatically during the loading phase. During ventral furrow formation, the leading edge of the ectoderm travels toward the midline at a rate of approximately 0.1 µm/s (He et al., 2014). Since this is slower than the deformation rate in our experiments, dissipative forces such as friction must also be lower during ventral furrow formation. Therefore, we can conclude that tissue deformation is also largely adiabatic during this morphogenetic process.

The time course of deformations and forces observed in our pulling experiments also reveal additional key mechanical features of the cellularizing epithelium.

One important feature is that the tissue is significantly elastic. During the loading phase (in which constant displacement is imposed), the applied force increases linearly with displacement, rather than with the rate of displacement (Figure 3c). If viscous forces were dominant then the force should be proportional to the rate of displacement, which is not the case. During the unloading phase, after the cantilever has been withdrawn, the tissue relaxes largely, though not entirely, back to its initial configuration (Figure 2a). These observations are both clearly consistent with the tissue itself being either an elastic or viscoelastic material.

The arguments given in the previous section for adiabatic deformation also mean that elastic forces dominate over friction and viscous shear forces in this system during the loading phase. If friction and viscous shear forces contribute only in a minor way during loading, what is primarily limiting tissue deformation in response to applied force? The observation that force on the cantilever drops by only about half (not to zero) during the holding phase suggests that boundary constraints may be important. To see why, consider a linear elastic spring which is immersed in a highly viscous fluid, and which is free at both ends. In response to force applied at one end, the trailing end will also move, though more slowly due to drag, and the spring will stretch. If the leading end is then held in a constant position, the trailing end will eventually catch up, the spring will eventually reach its resting length, and thus the spring will eventually cease applying a backward force on the leading end. If, however, this spring is fixed in place at its trailing end, a similar sequence of pulling then holding the free end will yield a different force profile. The spring will stretch while the free end is being pulled, but when the free end is then held fixed, the trailing end is unable to catch up. The spring should eventually reach a state of uniform deformation, but it cannot reach its resting length, and so it will continue to exert a force on whatever is holding the leading end in place. This analogy suggested to us that the reason the force on the cantilever never approaches zero during the holding phase of our pulling experiments is due to some type of boundary constraint. The nature of this ‘boundary constraint’ was not obvious to us at first, but will be explored in the computational studies below.

With these mechanistic insights in mind, we wanted to build a computational model both to test our preliminary conclusions and to estimate the elastic modulus of the cellular layer. To represent our experimental system computationally, we developed a three-dimensional (3D) model representing the entire embryo (Figure 4), in which geometries of all the relevant structures were taken from literature (Mavrakis et al., 2008) and measurements (Figure 4b). The cells are represented as open prisms with lateral and apical surfaces modeled as networks of elastic springs (each with spring constant k); the basal surfaces are left open (Figure 4a''). The cytoplasm filling the cell interiors and the yolk are modeled as a Newtonian fluid, with viscosity set to the experimentally measured value of 1000 cP (Doubrovinski et al., 2017; Selvaggi et al., 2018). The surface of the embryo is a rigid shell, with a no-slip condition imposed, representing the vitelline membrane, and the perivitelline fluid between the vitelline membrane and the cellular layer is approximated by a Newtonian fluid with the viscosity of water (~1 cP) (Figure 4b). External force is applied to a single cell edge, representing the effect of the cantilever.

Figure 4 with 2 supplements see all

Download asset Open asset

Although the elastic components as described above are the springs, there is no obvious correspondence between individual springs and any particular cellular structure. In other words, a spring constant derived from this model would not reflect the mechanical properties of any particular cellular component. However, there is biological significance to the triangular network of springs which comprise a single cell surface; it behaves like an elastic shell or sheet, and is a good model for the biological cell surface (comprised of plasma membrane and its associated cytoskeletal network). Previous work (Seung and Nelson, 1988) has shown that in an elastic sheet modeled as a triangular network of springs, the elastic (or Young’s) modulus of the sheet (E) and the spring constant (k) obey the following relationship: $E=2k/\sqrt{3}$ . Therefore, elasticity will be represented by either E or k in the following discussion, depending on the context. Note that for this type of approximation, the Poisson’s ratio is 1/3.

Technical details on the numerical implementation of our model are deferred to the Appendix. However, we briefly summarize key aspects of our modeling approach here. We used the immersed boundary method (Peskin, 2002; Peskin and Printz, 1993), which describes (1) how fluid responds to forces applied on it by an immersed solid object and (2) how the solid moves in response to external forces and fluid-solid interactions.

For the fluid domains (in this case, the cytoplasm and perivitelline fluid), we assume that the flow of the fluid is governed by Stokes equations. As a result, velocity of the fluid at any spatial location is uniquely determined from external forces driving the flow. In practice, to calculate this velocity distribution, we discretize the fluid domain by subdividing it into small tetrahedra with fluid velocities specified at every vertex (or node). The Stokes equations provide a relation between external forces and fluid velocities (Happel and Brenner, 1973). Thus, if we specify three spatial components of external force at N fluid nodes, we can obtain three components of velocities at each one of those N nodes by solving the discretized version of the Stokes equations.

The external force on the fluid comes entirely from the movement of the immersed solids, which this case consists of cell membranes. In this particular model, we need to calculate how the dynamics of the cell membranes are driven by the external force applied by the cantilever. As explained above, we model each cell membrane as a triangular network of springs. This network is in turn surrounded by fluid nodes, describing the fluid.

Each simulation step consists of five ordered sub-steps. (1) First, one calculates forces at the nodes of the triangulated elastic network. For most nodes, this is simply the sum of forces from elastic springs adjacent to that solid node; however, for nodes in contact with the cantilever, the force from the cantilever is also included. (2) Next, forces from the solid nodes are transferred to the fluid. The force on each fluid node is calculated from the force of the nearest solid nodes, with a correction for distance (see Appendix for details). (3) Now that external force on the fluid (due only to mechanical interactions with the immersed solid) is specified, velocities at fluid nodes may be calculated from the Stokes equations. (4) Next, velocities are transferred from fluid to solid. Specifically, the velocity on each solid node is calculated from the velocities of the nearest fluid nodes, with a correction for distance. (5) Next, coordinates of all solid nodes are updated by translating solid nodes with their now known velocities for a small interval of time dt. This modifies the configuration of the elastic solid and forces at all solid nodes must be recalculated anew. Iterating these steps, time evolution of the system is determined for some desired interval of time.

During the unloading or relaxation phase, the computational model described above has only one free parameter: the spring constant, k, of the individual springs (or equivalently, the Young’s modulus, E, of the cell surfaces). We could thus run the model with a range of values for k and test which values, if any, fit the experimental data (Figure 5). To do this, we fixed a value for k and imposed a rate of tissue deformation at a single edge which matched that imposed experimentally during the loading and holding phases. We then removed the applied force and allowed the model tissue to relax (Figure 5a). To quantify the relaxation rate, we calculated an exponential decay rate τ from these data. Carrying out this parameter sweep, we found that a value of k=0.001 nN/μm (or equivalently E=0.0012) resulted in a curve that matched the data extremely well (Figure 5a and c, Simulation 1).

Figure 5

Download asset Open asset

Having found that a spring constant of 0.001 nN/μm allowed us to realistically model experimental data from the unloading (or relaxation) phase of our experiments, we next wanted to test whether our model could also accurately capture the loading (or pulling) phase. During the loading phase, we had imposed a predetermined rate of displacement for a single edge, so by design the displacement matches that seen in vivo. However, in the loading phase we have an additional quantity that comes into play: the applied force needed to drive the tissue deformation. The approach we used to determine the instantaneous force needed in our model to generate the predetermined displacement is described in the Appendix. For each value of k examined in our parameter sweep, we determined the force required for the tissue deformation and plotted this predicted force as function of time (Figure 5b). The predicted force increases linearly during the loading phase, which agrees well with the experimentally observed linear increase in force during loading (Figure 3c). However, for the value of k which best matched the decay dynamics (k=0.001 nN/μm, Simulation 1), the predicted maximal force at the end of the loading phase was less than 0.2 nN (Figure 5b and c, Simulation 1). This value is an order of magnitude lower than the forces measured directly from the deflection of our calibrated cantilevers at the end of either the loading (~11 nN) or holding (~5 nN) phase. Conversely, for the value of k (0.050 nN/μm, Figure 5, Simulation 3) for which the predicted maximal force (3.4 nN) most closely resembled the experimental data, the decay dynamics were far too rapid, with τ=8.9 s, compared to τ=60 s in the experimental data.

What is the reason for this discrepancy? Our predicted loading force was based on the assumption that the dominant effects on tissue mechanics (as captured in our model) are the same in the loading and unloading phase. Therefore, we began to consider ways in which this assumption might not be true. We thus began a series of computational experiments using simplified models to explore how specific features of the model affect loading and unloading.

One feature we examined was the effect of network topology on the dynamics. Recent work has explored the physical consequences of the so-called ‘floppy’ networks, in which the number of degrees of freedom exceeds the number of constraints, versus ‘stiff’ networks, in which constraints exceed degrees of freedom (Broedersz and MacKintosh, 2014; Wyart et al., 2008). For simplicity (and to reduce simulation time), we used a model tissue in which the epithelial layer was replaced by a simple network of elastic springs, and viscous forces from the cytoplasm and yolk were replaced by a constant drag force (see Materials and methods). We then explored the different effects of stiff versus floppy networks by using either a triangular (stiff) or hexagonal (floppy) network of springs.

During the initial loading phase, we imposed a fixed deformation rate at a chosen point, and determined the instantaneous external force required to achieve this deformation. We performed this simulation for models incorporating either triangular or hexagonal networks, and also varied the spring constant k assigned to the springs or edges. For a given spring constant, the force required to deform the triangular network was greater than the force required to deform the hexagonal network by the same amount (Figure 6a). In spite of this quantitative difference, these networks responded to pulling in a similar way qualitatively: in both cases, the force needed for a given deformation increased as the elastic modulus increased (Figure 6a).

Figure 6

Download asset Open asset

In the subsequent unloading phase, we removed the externally applied force and observed how the tissue relaxed over time (Figure 6b). In the triangular cases, the tissue relaxed back to its initial configuration more rapidly as k increased (Figure 6b and c). This agrees with the intuition that tissue with higher elasticity should spring back more rapidly. In the hexagonal-mesh case, however, the tissue relaxed back to its initial position more slowly as k increased (Figure 6b). Furthermore, the extent of tissue relaxation was also lower (Figure 6b and d). These results indicate that the response of hexagonal and floppy networks during unloading is qualitatively different.

There is an obvious parallel between these types of networks and different compartments of the cellularization-stage epithelia in our experiments. The apical surface appears most similar to a stiff elastic sheet, while the lateral surfaces and basal edges form a semi-hexagonal, and therefore floppy, network. Since the apical and basal-lateral compartments would therefore be expected to contribute differently to the loading versus unloading phases of our experiments, we hypothesized that our earlier mismatch of maximal forces calculated from the loading versus unloading phase may be due to a difference in the elastic modulus of the apical and basal-lateral compartments.

To test whether different compartments of the cellularization-stage epithelia exhibit different elastic properties, we modified an experimental setup we had used previously (Doubrovinski et al., 2017). In those previous experiments, we injected a 20–30 µm diameter ferrofluid droplet into the embryo, positioned it in the cellular layer, and used a magnet to apply force across the basal-apical extent of the cell. In these previous experiments (reproduced in Figure 7a–a'), the speed of deformation decreased as the droplet approached the magnet, which we attributed to a build-up of elastic stress in the epithelium. In our new experiments, we used a smaller (approximately 10 µm diameter) ferrofluid droplet, which we could then position entirely inside the cell, allowing us to specifically apply force to a small region at the intersection of the apical and lateral surfaces of the cell, without directly contacting the basal edges (Figure 7b). In these experiments, the speed of deformation was constant or even increasing as the droplet approached the magnet, demonstrating that elastic stresses in this case were very small compared to the pulling force (Figure 7b'). This difference in dynamics demonstrates that the apical and lateral portions of the cell are softer, that is have a lower elastic modulus, than the basal edges.

Figure 7 with 1 supplement see all

Download asset Open asset

We then returned to our full computational model of the embryo (Figure 4), this time assigning different values of the spring constant k in different cellular compartments. To test the effect of varying elasticity in the floppy cellular components, we did a series of simulations in which we held k_apical fixed and assigned increasing values to k_lateral and k_basal (Figure 7c–c''). Increasing k_lateral and k_basal 20- to 30-fold resulted in much higher values for maximal force F_max (going from 0.2 nN to nearly 7 nN), but resulted in very similar values for the decay constant τ (~50 s). This validated our hypothesis that increased elasticity within floppy structural components should primarily affect mechanics during loading (characterized by F_max), with only minor effects on unloading (characterized by τ). In contrast, when k_basal and k_lateral were kept approximately constant, but k_apical was increased about sixfold, this resulted in a dramatic change of decay constant τ, from 39 s to 9 s (Figure 7—figure supplement 1).

Notably, one of the parameter sets used in our simulations (k_apical = 0.007, k_lateral = 0.126, k_basal = 6.825) resulted in a maximal force (7 nN) that fell within the 5 nN to 11 nN range of maximal force measured directly from cantilever deflection in experiments. Furthermore, despite the large values for basal/lateral elasticity (compared to k=0.05 in Figure 5), the fit to relaxation dynamics remained relatively good. To summarize, assuming that basal domains are much stiffer than the apical and the lateral domains suffices to reconcile the model with our measurement data. Note that this key assumption is independently verified in our ferrofluid experiments (Figure 7a and b) and that their interpretation in no way relies on either the cantilever-based experiments or their computational analysis. Note additionally that the final estimated values should be considered order-of-magnitude estimates, given the present measurement accuracy.

It is important to note that we expect that different distributions of elasticity among the floppy portions of the cell would yield similar results. Since we currently do not have experimental tools to measure, and thus computationally constrain, the ratio of lateral-to-basal elasticity, we cannot at this point precisely determine those numbers independently. However, our computational and experimental results combined do lead us to two firm conclusions. (1) The elastic modulus we estimated from relaxation dynamics (E approximately = 0.001–0.01) provides a good estimate specifically of apical elasticity. (2) Basal elasticity in this tissue is substantially higher than apical elasticity.

By examining our computational model, which closely reproduces the dynamics seen in vivo, we found ourselves able to answer a question that had puzzled us earlier. Why do tissue dynamics appear to be subject to boundary effects if the embryo is an ellipsoid and, thus, the cell layer appears not to have boundaries? One possibility is that some feature of the epithelium (.e.g., the topology of the cell network) stabilizes the tissue structure in a way that cells some distance from the cantilever would be immobile, forming a sort of pseudo-boundary. Plotting the displacement of cells in our model reveals that this explanation is not true. Even at short pulling times, cells throughout the model embryo move from their initial positions, even on the opposite side of the embryo from the force application, and even at the distant anterior and posterior poles (Figure 8a–a', Video 3). We next examined previously acquired experimental data (Doubrovinski et al., 2017) to see if the same might be true in vivo. Indeed, when concentrated force was applied to one side of the embryo using a ferrofluid droplet, cells on the opposite side of the embryo clearly moved as well (Figure 8b–b', Video 4). Although cells appear to be moving in a coherent direction (e.g., counterclockwise clockwise in Figure 8b–b') along the sagittal plane, the tissue is clearly not freely rotating, since the displacement drops off away from the applied force. A consideration of embryo shape suggests two obvious reasons why it does not rotate freely under applied force. Cells moving between the highly curved poles and the flatter lateral regions would need to change shape and/or neighbor relationships, and there may energetic barriers to this transition. In fact, our simulations reveal an upregulation of stress near the poles (Figure 8c). Additionally, if cells cannot rearrange completely freely, an imaginary ring composed of cells encircling the embryo crosswise must deform into a somewhat larger ellipse under the described tissue motion (Figure 8d), which would be expected to result in increased elastic tension among these cells. In summary, geometrical constraints preventing free rotation of the tissue allow for the build-up of elastic stress throughout the tissue when force is applied locally.

Figure 8

Download asset Open asset

Video 3

Download asset

Video 4

Download asset

All experiments used the membrane-Cherry fly line described in Martin et al., 2010. The fluorescent signal proved too faint for our purposes, so a membrane dye was used to visualize membrane outlines instead (see below).

Embryos were collected using on grape juice agar (Genesee Scientific) plates with a smear of yeast paste. To select embryos of a suitable stage, plates were covered in light halocarbon oil (Halocarbon Oil 27, Sigma). Embryos with a distinctive faint halo in the periphery, indicative of early cellularization, were transferred to small square mats of lab paper towel for further manipulations. The embryos were dechorionated in bleach, then quickly rinsed 10–12 times in distilled water. To attach the embryos to a slide, a thin streak of embryo glue (a saturated solution of adhesive from double-sided Scotch tape dissolved in heptane) was placed across the length of a coverslip and allowed to dry for 5–10 min. Dechorionated embryos were picked up and placed on top of the glue streak using an eyelash. A total of 8–12 embryos were chosen for each round of imaging. The embryos were then desiccated in a closed container with Drierite (Fisher Scientific) for 22 min before further manipulations.

Needles used for dye microinjections were fabricated from borosilicate capillaries (B100-75-10, Sutter Instruments) using a pipette puller (P-2000, Sutter Instruments). Pipettes were backfilled with 25-fold diluted CellMask Deep Red Plasma Membrane stain (Thermo Fisher Scientific) using a gel loading tip (Eppendorf). Pipettes were broken at the tip by gently nudging the pointed end against an edge of a coverslip, under a drop of oil.

To inject dye, the tip of a pipette was carefully pushed against the vitelline membrane so as to insert it into the perivitelline space without piercing the cellular layer. An Eppendorf injector (FemtoJet 5247, Eppendorf) was used with the following settings: compensation pressure: 10 hPa, injection pressure: 220 hPa. Typically, the injection was performed using the ‘clean’ command such that maximal pressure was applied.

After injection, the embryos were cut along the side opposite to the site of injection for subsequent cantilever insertion.

A glass pipette was pulled the same way as done when making microinjection needles. The pipette was then similarly fractured at the tip by gentle contact with a glass coverslip, but in air with no added oil. To make the polymer filling for the glass needle, PDMS (Sylgard 184, Dow Corning) was mixed with a stock solution of BODIPY (BODIPY 493/503, Thermo Fisher Scientific; 1 mg per ml stock solution in anhydrous DMSO) in proportions 3:2 (by volume). The mixture was placed in a dish of aluminum foil and left overnight at 80°C to evaporate off DMSO. Next, the PDMS-BODIPY mixture was mixed with PDMS crosslinker in proportions 10:1 by weight. This mixture cures within 48 hr if left at room temperature, but may be stored for months at –20°C without curing. A small droplet of PDMS-BODIPY-crosslinker mixture was placed in the center of a weighting boat. A glass pipette is placed into the weighting boat so as to immerse its (open) tip in the droplet, simultaneously leaning the other end of the pipette against the edge of the weighting boat. The weighting boat is left at 80°C for 5 hr allowing PDMS to cure. Finally, the pipette is carefully lifted out from a (now cured) PDMS droplet. In this way, the inner part of the pipette tip is filled with crosslinked fluorescent PDMS.

The PDMS-filled pipette was bevelled using a micropipette grinder (E-45, Narishige) for 20 s with the speed dial set to 7. This is needed to remove the very end of the tip since otherwise the very tip of the cantilever will be too thin and will tend to curl on itself.

Finally, the outer layer of glass was etched away using Armour Etch (Armour Glass Etching Cream, 2.8 Oz), a paste comprising a mixture of hydrofluoric acid and fluorites. A spectrophotometer cuvette was filled with Armour Etch cream halfway, and a layer of light silicone oil (Silicone oil for oil baths, Sigma) was added on top to prevent solvent evaporation. The tip of the pipette was immersed into the etching paste and was left to etch for 3–5 min. To observe the pipette tip through a vertically mounted microscope, a retrofitted MicroForge (MF-900) with a custom mount for the cuvette was used (Figure 1—figure supplement 1d). This completely dissolved the exposed glass layer to expose the core of cured PDMS. The pipette was then lifted out from the cuvette and gently washed with a solution of Triton X-100 and saline (0.4% NaCl and 0.03% TX-100 in distilled water). Finally, the pipette is bent 90 degrees in flame of a cigarette lighter to simplify its mounting in a pipette holder.

It is important to note that filling capillaries from the front, by capillary action, is absolutely critical for generating cantilevers of reproducible properties. When microfabricating cantilevers using pulled glass pipettes as a mold, one may try introducing PDMS into the pipette by backfilling (as we in fact did initially). When microcantilevers are fabricated in this way, their Young’s modulus (measured by subjecting the probe to viscous drag) varied widely. To track the origin of this undesired variability, we performed focused ion beam (FIB) milling on our microfabricated probes. In the course of milling, half of the cantilever beam was milled away along its length so as to expose and examine the mid-cross-section of a sample using scanning electron microscopy. It was found that portions of the cantilever were in effect hollow, due to sub-micron-sized bubbles trapped in the polymer (Figure 1—figure supplement 1a). We believe that bubbles were formed due to the air trapped inside the capillary tip when it is being backfilled. In contrast, when capillaries are filled from the front using capillary action, FIB milling reveals no bubbles trapped within the polymer (Figure 1—figure supplement 1b), and measurements of the Young’s modulus are consistent.

The glass pipette with its PDMS cantilever tip was translated using a piezo actuator (PXY200SG two-axis piezo actuator, Newport). The piezo was controlled with an amplifier (NPC3SG three-channel piezo amplifier, strain-gauge position control, Newport) interfaced with MATLAB. The MATLAB code used to control the piezo is available for downloaded from https://github.com/doubrovinskilab/cantilever_embryo_rheology (copy archived at Cheikh, 2023). To mount the piezo actuator on a microscope stage, an optical breadboard was screwed to the stage using two custom threaded holes. A three-axis stage (PT3 - 1" XYZ Translation Stage with Standard Micrometers, Thorlabs) was mounted on the breadboard. The piezo actuator was attached to the three-axis manual stage using two 45-degree angle brackets and a custom adopter mount.

Confocal laser scanning microscopy was done using a Zeiss LSM 700 with excitation wavelengths of 488 nm (green) and 639 nm (far red).

The computational methods used in this paper (including cantilever calibration, models, and numerical implementation) are described in detail in the Appendix.

In order to obtain the Young’s modulus of a cantilever, the cantilever was dragged at a constant speed of 53 µm/s through Tween 20, a fluid whose density is density ${\displaystyle \rho ={10}^3\mathrm{k}\mathrm{g}/{\mathrm{m}}^3}$ and viscosity ${\mu }_f=460\mathrm{c}\mathrm{P}$ (information provided by Sigma-Aldrich). Using these data and the measured dimensions of the cantilevers calibrated, the characteristic Reynolds number can be estimated to $Re~{10}^{-5}$ , implying that inertial effects are negligible. The dynamics of the cantilever is thus governed by the balance between viscous and elastic forces. Viscous drag was calculated based on the slender body theory following Cox, 1970, which is justified by the small diameter-to-length ratio of the cantilever. In particular, the drag force per unit length reads:

where $R\left(x\right)=R_{min}+x(R_{max}-R_{min})/L$ and $R_a={(R}_{max}+R_{min})/2$ are the local and average radii of the tapered cantilever, respectively, $h$ is vertical deflection at a given position ($x$) and time ($t$), and $L$ is the length of the cantilever. On the other hand, since the cantilever length remains constant, we use inextensible beam theory to obtain the expression for the force density due to bending (Semler et al., 1994):

where $I\left(x\right)$ is the area moment of inertia (at position $x$). To determine the value of the Young’s modulus $E$, we performed a parameter sweep varying the value of $E$. To do this, we solve the equations following from the force balance (along the direction perpendicular to the axis of the cantilever and assuming small displacement): $F_v+F_e=0$. This was done using the open solver FreeFEM (Hecht, 2012).

Here, we describe the numerical method used to simulate tissue deformation. Our approach is a combination of the immersed boundary method (IBM) (Cuvelier et al., 1986; Peskin, 1972; Peskin and Printz, 1993) and the finite element method (FEM) (Boffi et al., 2007). The chosen methodology is motivated by the requirement to describe fluid-structure interactions when treating the dynamics of embryonic cells that are elastic on the surface and filled with viscous fluid in the interior (Doubrovinski et al., 2017).

The IBM approach involves considering an elastic solid immersed in a fluid as shown in Figure 4—figure supplement 1. In what follows, ${\Omega }_f$ denotes the region occupied by the fluid with boundary ${\Gamma }_f$ , while ${\Gamma }_s$ is the interface of the fluid-immersed elastic solid. Initially, the fluid domain ${\Omega }_f$ is discretized into finite elements (Figure 4—figure supplement 1b-d). Within the IBM approach, the geometry of this (Eulerian) grid remains fixed throughout a simulation. The interface of the solid is meshed into triangles (whose edges will be modeled as linear springs, see Figure 4—figure supplement 1b). The nodes of the corresponding discretized (Lagrangian) solid will be translated according to the relevant equations of motion which we turn to next.

Describing the evolution of the solid interface (the cell surfaces) requires simulating the ambient fluid, whose dynamics are given by the Stokes equations:

Here, µ is dynamic viscosity, while $u_f\left(x,t\right)$ and $p_f\left(x,t\right)$ are the velocity and the pressure fields, respectively. Force density acting on the fluid is $f_f$ . In our case, $f_f$ derives entirely from forces originating from the immersed structure ${\Gamma }_s$ . We are justified in neglecting inertia since a typical Reynolds number in our problem may be estimated to be on the order of ${\displaystyle Re={10}^{-8}}$ . More specifically, in a typical experiment, embryonic tissue is subjected to concentrated force as the tissue is pulled by a cantilever. The typical speed $U$ with which the cantilever is translated was $0.005\mathrm{}\mu m/ms$. The typical size (length) of an embryo $L$ is $500\mathrm{}\mu m$. The viscosity of the embryonic cytoplasm was previously measured to be on the order of $1.0\mathrm{}nN\bullet ms/{\mathrm{\mu }m}^2$ (Doubrovinski et al., 2017; Selvaggi et al., 2018). The density of the cytoplasm is on the order of that of water $\rho ={10}^{-6}\mathrm{\mu }g/{\mathrm{\mu }m}^3$ . Hence, we obtain ${\displaystyle Re=\frac{\rho UL}{\mu }=2.5\times {10}^{-8}}$ as above.

To describe the surfaces of model cells, we assume that the cortex may be described as a thin linearly elastic shell characterized by a (2D) Young’s modulus and a Poisson’s ratio.

Finally, the dynamics of the solid are coupled to that of the fluid by imposing a no-slip boundary condition, such that the velocity of the fluid matches that of the solid at the solid interface. The net force density on the solid ($f_s$ , which is defined below) is equated with the force on the fluid (through the $f_f$ term). Note that these are the key standard assumptions exploited in a typical IBM implementation (Peskin, 1972; Peskin and Printz, 1993).

For the fluid, we start by casting the governing equations for the fluid into their weak form:

where $\psi$ and $\varphi$ are test functions for the pressure and the velocity field, respectively (Boffi et al., 2007; Boffi and Gastaldi, 2003; Cuvelier et al., 1986). The spatial domain is discretized into P1b-P1 finite elements which satisfy the inf-sup condition (Boffi and Gastaldi, 2003; Boffi et al., 2007; Cuvelier et al., 1986). More specifically, we choose the cubic P1b element for the velocity test function ${\displaystyle \varphi }$ and the linear P1 element for the pressure test function $\psi$. Note that there are $N$ discrete nodes for the velocity field and $M$ nodes for the pressure field in what follows.

Using the above discretization, the relevant fields are approximated as:

In terms of the discretized variables the Stokes equations yield the Galerkin form:

Equations 6 and 7 define a finite element system of 3N+M equations and 3N+M unknowns (in the 3D case). Assuming that the fluid boundary is Dirichlet, the linear system will take the standard form

where

The plasma membrane surfaces of the embryo can be thought of as a set of interconnected linearly elastic 2D plates. A particularly time-efficient implementation for representing elastic plates involves discretizing each plate as a set of interconnected triangles whose edges are linear springs (Seung and Nelson, 1988). More specifically, it may be shown that in the limit of small deformation, the discrete triangulated mesh approximates a 3D linearly elastic material whose Young’s modulus is given by $E=2k/\surd 3$ (with $k$ being the spring constant of the mesh) and Poisson’s ratio ${\displaystyle \sigma =1/3}$ (Seung and Nelson, 1988). Therefore, our immersed solid ${\Gamma }_s$ is modeled as a triangular mesh of springs. The no-slip boundary condition at the fluid-solid interface amounts to moving the solid nodes with the local velocity of the fluid:

where ${\displaystyle X_s}$ is a given point on the surface of the solid. Every spring comprising the discretized solid contributes a force dipole at its two adjacent nodes:

where $F_s$ is the magnitude of the elastic force at a solid node, $k$ is the spring constant, ${\displaystyle \mathrm{\Delta }L}$ is the difference of the corresponding instantaneous spring length and its rest length, while $\widehat{s}$ is a unit vector along the respective spring.

The coupling of solid and fluid dynamics is implemented within the IBM framework in a standard way (Peskin, 1972; Peskin and Printz, 1993). In brief, force density of the solid is transferred to the adjacent fluid nodes (the ‘spreading’ step), contributing the $f_f$ -term in Equation 2. Next, the velocity and pressure fields of the fluid are obtained from the system of linear equations (Equation 8). Finally, in the ‘interpolation’ step, the velocity of the fluid is interpolated to the solid nodes and those nodes are then advanced through a time-step using a Euler-forward scheme.

In this section, we validate our numerical solver by comparing simulation results to the solution of an analytically tractable fluid-structure interaction problem. The problem, illustrated in Figure 4—figure supplement 2, involves stretching an elastic plate submerged in the middle of a viscous fluid. The plate is fixed on the left, with the right end subjected to a constant force.

Neglecting inertia, we can write the equations governing the stretching of the plate as:

where $\delta \left(x,t\right)$ is the displacement of the plate, $E$ the 2D Young’s modulus of elasticity, µ the viscosity of the fluid, and $\frac{\partial \delta \left(x,t\right)}{\partial \mathrm{t}}=u(x,H,t)$ is the $x$-velocity component of the fluid. Equations 13 and 14 are the boundary conditions for the problem, Equations 15 and 16 are the initial conditions, and Equation 12 represents the balance of forces acting on the plate (equating stress in the fluid to force density on the solid). The factor of 2 in Equation 12 comes from the fact that viscous drag is present on both sides of the elastic plate.

To solve this system analytically, we assume the distance between the plate and fluid boundaries is small in comparison to the length of the plate (i.e., $H\ll L$). In this case, the flow is well approximated by the Poiseuille solution such that the $y$-velocity is identically zero, whereas $x$-velocity $u(x,H,t)$ is linear in between the elastic plate and the horizontal walls of the fluid domain. With these assumptions, Equation 12 may be recast as a diffusion problem

where $\eta =2\mu /H$ and $D=E/\eta$ is the diffusion coefficient. The solution of this problem is given by the Fourier expansion

which may be considered a sum of the steady-state solution and a transient. This expression, which is the analytical solution to the problem, is shown as the red line in Figure 4—figure supplement 2b-c.

We also simulated the problem using our numerical scheme to obtain an approximate solution. This simulated solution is shown as the blue line in Figure 4—figure supplement 2b-c.

The analytical expression and numerical simulation are in very close agreement. Moreover, as expected, the agreement systematically improves with decreasing thickness of the fluid layer (data not shown). Taken together, simulation results in Figure 4—figure supplement 2 provide a rigorous validation of the chosen numerical scheme.

In this section we detail simulations used to interpret the pulling experiments probing the rheology of embryonic epithelium. The geometry of a model Drosophila embryo is illustrated in Figure 4. The computational domain is an ellipsoid whose dimensions approximately match the well-documented dimensions of an embryo (Mavrakis et al., 2008). The outer solid boundary represents the rigid vitelline membrane. Beneath the vitelline membrane is a thin (3 μm) gap of perivitelline fluid. The rest of embryonic interior is described as a uniform Newtonian fluid with a viscosity of 1000 cP, as previously measured (Doubrovinski et al., 2017; Selvaggi et al., 2018). Cells are hexagonal prisms whose faces are elastic plates modeled by a set of inter-connected springs (Figure 4a–a''). Dimensions of cells are based on experimental data.

A P1b-P1 FEM mesh was generated so as to discretize the whole computational domain (including the model perivitelline space, cytoplasm, and yolk). The mesh was further refined in the region corresponding to the perivitelline fluid. Dirichelet (no-slip) condition was prescribed at the outer domain boundary corresponding to the vitelline membrane.

To capture the cantilever pulling experiments, simulation is divided into three phases:

1. In the loading phase, one cell edge is subjected to force a pulling force for 60 s. In our experiments, the imposed speed of the cell edge that was being pulled remained approximately constant during the course of loading. To match the measurements, in our simulations, the external pulling force was applied through a proportional control system (P-controller) such that the externally applied force is given by

   • (19) ${\overrightarrow{F}}_{external}\left(t\right)=K_p\overrightarrow{e}\left(t\right)-{\overrightarrow{F}}_{internal}\left(t\right)$

   • $K_p$ was set to be 0.2 nN/μm; ${\overrightarrow{F}}_{internal}$ represents the sum of the internal forces on the cell edge being pulled; the error value $\overrightarrow{e}\left(t\right)$ is the difference between the measured and the simulated displacement of the cell edge that is being pulled. The external force (e.g., Figure 7c') in our simulations was obtained from Equation 19.

2. In the holding phase, the cell edge was held fixed for an additional time of 60 s.

3. In the unloading phase, the external pulling force was removed.

All simulation code is publicly available under https://github.com/doubrovinskilab/cantilever_embryo_rheology (copy archived at Cheikh, 2023).

To better understand simulation results obtained from the model in the previous section, we considered the behavior of floppy and stiff networks in the experimentally relevant geometry (Figure 6). The full model is extremely computationally intensive; for example, Simulation 4 in Figure 7c' required 3 weeks of simulation time using 200 processor cores. Therefore, we used a simplified model to allow us to more quickly explore parameter space. One simplification was including 2D spring networks confined exclusively to the surface of the ellipsoid representing the surface of the embryo (using a soft force of constraint). The more significant simplification in terms of computational time was that we replaced explicitly modeled Stokes fluid with effective drag.

As discussed below, the results from the simplified treatment cannot be quantitatively compared to the results from the full treatment including hydrodynamic interactions. However, to only gain qualitative insights into the behavior of our system, it is admissible to replace Stokes fluid dynamics with viscous drag, especially considering that we are mainly interested in the behavior of the solid and not the fluid. In our simulations (Figure 6), the networks were confined exclusively to the surface of the ellipsoid representing the surface of the embryo using a soft force of constraint. We checked that increasing the magnitude of the force of constraint by an order of magnitude did not affect the dynamics. All other aspects of these simulations were the same as in the case of the full hydrodynamic model.

For an example of how viscous drag can provide an inaccurate approximate for the physics modeled by Stokes fluid dynamics, consider the following example. A thin 2D disc moving through a 3D ambient Stokes fluid experiences force that is proportional to the radius; in contrast, using standard assumptions for drag, the force would be proportional to the area of the disc (Happel and Brenner, 1973). The use of drag in this example only approximates the results obtained from explicitly modeling the fluid in the limit where the fluid on either side of the disc is thin and bounded by another solid surface nearby (Happel and Brenner, 1973). Since the embryo is filled with a fully 3D fluid, we therefore chose the more physically correct model for our main simulations in spite of the limits this placed on our ability to efficiently explore large regions of parameter space.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Jeffrey Woodruff for helpful discussions. We also thank Salena Fessehaye for discussion of the manuscript. We gratefully acknowledge the help of the Characterization Center for Materials and Biology at UT Arlington, in particular Dr. Jiechao Jiang and Prof. Efstathios I Meletis, for their help with AFM characterization of the probe. We thank the Nano3 facility at UCSD for their help with the cryo-EM aspects of the project. We thank Saiana Khandarkhaeva, Natalia Dubrovinskaia, and Leonid Dubrovinsky for their help with the FIB modification of the commercial AFM probes. The work was supported by the Robert A Welch Foundation (grant I-1950-20180324) and the NIH (grants 1R01GM134207 and 1R21HD105189).

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.