{'fraction_non_alphanumeric': 0.0334633220910624, 'fraction_numerical': 0.01046058178752108, 'mean_word_length': 4.71839686605394, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Encoding Mechano-Memories in Actin Networks

Understanding the response of complex materials to external force is central to fields ranging from materials science to biology. Here, we describe a novel type of mechanical adaptation in cross-linked networks of F-actin, a ubuiquitous protein found in eukaryotic cells. We show that shear stress changes its nonlinear mechanical response even long after that stress is removed. The duration, magnitude and direction of forcing history all impact changes in mechanical response. The `memory' of the forcing history is long-lived, but can be erased by force application in the opposite direction. We further show that the observed mechanical adaptation is consistent with stress-dependent changes in the nematic order of the constituent filaments. Thus, this mechano-memory is a type of nonlinear hysteretic response in which an applied,"training"strain modifies the nonlinear elasticity. This demonstrates that F-actin networks can encode analog read-write mechano-memories, which can be used for adaptation to mechanical stimuli.

Biological cells sense mechanical stimuli and use these cues to control their physiology 1 . Elucidating such mechanosensory mechanisms in living matter has implications for both cell biology and smart materials design 2 . While molecular-scale mechanotransduction mechanisms are well appreciated 3,4 , how these are integrated to control mechanical response at cellular length scales is not well understood. Networks of the semi-flexible biopolymer actin are essential determinants of the mechanical behaviors of eukaryotic cells 5 . Previous work has shown that the mechanical response of biopolymer networks can be regulated either transiently 6 or irreversibly 7 by external force, necessitating mechanochemical feedbacks to enable long-lived and dynamic mechanical adaptation 8 .
When actin is polymerized in the presence of crosslinking proteins, space-spanning networks are formed (Fig.1a). The mechanical response of actin networks is viscoelastic 9-11 and highly nonlinear 12 . Typically, the strain stiffening of actin networks is reversible with no dependence on the history of applied stress 6 . However, polymeric systems of both biological and synthetic origin can also demonstrate irreversible work hardening 7 and other permanent rheological changes in response to applied stress [13][14][15] . The irreversibility of work hardening limits the potential for these materials to dynamically modify their mechanical properties or sense external force.
Here, we probe the rheology of a physiologically relevant network of actin filaments cross-linked with the protein filamin (FLN). FLN-crosslinked actin networks exhibit a dramatic strain stiffening [16][17][18] . FLN cross-links are also highly dynamic, with an off rate of ∼ 1 s −1 19 that results in a broad spectrum of stress relaxation 20 . We polymerize 24 µM G-actin in the presence of 5% FLN within a parallel plate rheometer and probe the mechanical properties via a symmetric, triangular strain 'read out' (RO) pulse ( Fig.1b). At low strain, a linear stress-strain relationship is observed but the stress increases nonlinearly at higher strains. The concavity in the stress-strain curve reflects the nonlinear strain stiffening; this behavior is symmetric for strains in the opposing direction due to the isotropic nature of the network (black triangles, Fig.1d). To explore how applied stress can alter this rheology, we apply a large constant stress in the positive direction (σ + = 4 Pa) for 250 s (Fig.1c). Upon removal of the stress, we wait for 50 s and read out the mechanical response by applying a RO pulse. Surprisingly, we observe a dramatic increase in the stress in the direction of the training with nominal reduction in stress in the opposite direction (Fig.1d, red squares); this indicates that the mechanical response is altered by the prior stress application. To explore the reversibility of this effect, we apply a constant stress of 4 Pa in the opposite direction (σ − ) for 250 s. Remarkably, after the negative training pulse, an enhanced stiffening is observed at negative strains (Fig.1d, blue circles) and all evidence of the previous positive training pulse is erased. To underscore the reversibility of this training, we apply successive cycles of positive and negative training. In Fig.1e, the stress measured between each cycle is plotted as a function of strain over six training cycles. These data demonstrate the ability to form and erase mechanical memories with robust reversibility. Such 'mechano-memory' is a type of nonlinear hysteretic response in which applied 'training' strain modifies the nonlinear elastic response of the material to subsequent loading.
The nonlinear mechanical response can be characterized by plotting the differential shear modulus, K = dσ dγ , obtained by differentiating the Lissajous curves numerically along the increasing strain magnitude, as a function of strain γ. For an untrained sample, K(γ) is symmetric around γ = 0 and shows stiffening observed at high  strains in both the positive and negative direction (black triangles, Fig.2a). For positively or negatively trained samples, the symmetry in K as a function of γ is clearly broken: K is enhanced for γ >0 for a positively trained sample (red squares in Fig.2a), whereas, for a negatively trained sample the enhancement happens for γ <0 (blue circles in Fig.2a). The parameter ∆K = K + −K − , where K + and K − are the moduli measured at strain magnitude 0.4 in the positive and negative directions, respectively, is used to quantify the strength of the memory formed (Methods). Due to the symmetry in the untrained sample, assignment of the direction of positive strain and stress is arbitrary. Once chosen, the sign convention is used consistently in all subsequent measurements for that sample.
To query whether the mechanomemory is long-lived, we vary the time (T w , Fig.1c) between removal of the training stress and the read out measurement. We find that, ∆K decreases weakly up to T w of nearly 1000 seconds ( Fig.2b), indicating that the training pulses encode a semi-permanent memory into the network. To examine how the training parameters can be used to modify the strength of the mechanomemory, we vary the magnitude and duration of the training stress. The duration of training tunes ∆K, with ∆K monotonically increasing from 2 to 8 Pa as training time T + is increased from 10 to 400 s (Fig.2c). Likewise, the training stress magnitude also tunes ∆K (see Fig.S1). For a sample maximally trained in the negative direction (σ − = 5 Pa, T − = 250 s), ∆K ≈ -8 Pa. When a positive training pulse σ + = 0.5 Pa is applied for T + = 100 s, ∆K increases to ≈ -5 Pa. When σ + > 4 Pa, ∆K > 0 (Fig.2d). After maximal training, it is possible to recover a mechanical response similar to an untrained sample with a suitable stress pulse in the opposing direction (Fig.S1).
To explore the effects of training stress, we study the time dependent changes in strain (i.e., creep) and recovery behavior. First, a training stress in the positive direction is applied for 300 s (Fig.3a). For a subsequent stress pulse in the same direction (σ + = 5 Pa, T + = 300 s), the resulting strain magnitude γ ++ increases to 0.8. After the stress is removed, the residual strain at 300 s, or γ r , is approximately 0 (Fig.3b). By contrast, when an identical stress pulse is applied in the direction opposite to the direction of training, the strain magnitude |γ +− | increases up to 1.2. After the stress is removed, γ r = 0.4 (Fig.3b). This increased residual strain reflects a plastic-like shape deformation upon stress reversal.
The differences in the creep observed in response to different directions of applied stress indicate that the system is effectively more elastic in the training direction compared to the opposite direction. Another way to think about this is that there is more energy dissipated in a creep measurement that occurs after reversal of stress direction. We directly calculate the dissipated power from the creep measurements and observe an an enhanced dissipation when the the direction of stress is changed (Fig.S2). The area under the power dissipation curves for these two different loading directions gives an estimate of additional energy that gets dissipated in reorganizing the system. We observe an increase in energy dissipation as a function of applied stress (Fig.S2). Due to the noise in calculating energy dissipation directly from creep measurements, we use an alternate strategy (Methods) of superimposing a small, oscillatory stress during the creep measurement to probe power dissipation (Fig.3c). Using this alternate strategy, we also measure an enhanced dissipation when the system is loaded in the direction opposite to that of training direction (Fig.3d). We then calculate the dissipative power required to rewrite a memory of similar strength but different sign (i.e., from +∆K to -∆K) and normalize it by that dissipated in the absence of a new memory formation (P +− /P ++ ) (Methods). We find that the power dissipation increases dramatically for stresses > 3 Pa and mostly contained within the first 100 s of stress application (Fig.3e). This underscores the stress-and time-dependent work done by the external mechanical stresses to dynamically alter the network mechanical response.
To elucidate the mechanism of mechanomemory, we compare the normalized differential modulus as a function of strain ( Fig.4a) for networks trained (σ + = 4 Pa) for a duration of T + = 0 s (untrained, black circles), T + = 200 s (upside down triangles) and T + = 400 s (squares) with a simple two-dimensional finite-element simulation (lines). In the simulation 2 µm long stiff filaments are cross linked by 100 nm long, flexible linkers to represent FLN (see Fig.S3). For the isotropic networks, the strain dependencies in the differential modulus arises primarily from non-linearities in the constitutive elements (e.g. cross-linkers and filaments) 21 .
Previous work has shown that small degree of filament alignment, or nematic order, leads to a dramatic change in the nonlinear elastic response 22,23 . Specifically, when a network of rigidly cross linked semiflexible filaments is sheared so that the compression direction lies along the nematic director, softening is observed at strains of less than 0.1% due to the buckling of the filaments under compression (Fig.4b, bottom schematic). When the same network is sheared in the opposite direction so that the same aligned filaments are now tensed, the system exhibits strain stiffening (Fig.4b, top schematic). These nonlinear elastic effects make the

Stress
FIG. 3. History dependent creep behavior and energy dissipation a, Schematic of the protocol to measure history dependent creep response: An untrained sample (ca = 24 µM actin + 5 % FLN) is first trained by σ+ = 5 Pa stress in the positive direction for 300 s. After waiting for 300 s, a stress of magnitude 5 Pa is again is applied for 300 s either in the training direction (positive) or opposite to it (negative) and then the stress is switched off (σ = 0 Pa) for 300 s. b, Temporal evolution of the magnitude of the creep strain when the stress pulse is applied in the training direction |γ++(t)| or opposite |γ+−(t)| to it. The residual strain (γr) estimated after 300 s of recovery is indicated for the case of strain evolution γ+−(t). c, Schematic of the protocol to estimate energy dissipation in rewriting a memory. d, An untrained sample (ca = 24 µM actin + 5 % FLN) is first trained by σ+ = 4 Pa stress in the positive direction for 250 s. Time dependent power dissipation is estimated (see Method) by superposing a small A.C. component of stress (amplitude 0.4 Pa) on the D.C. stress pulses (magnitude 4 Pa) applied in the same (P++(t)) or opposite direction (P+−(t)) w.r.t the initial training direction. The integrated area between P++(t) and P+−(t) indicated by the gray colored region gives the additional dissipated energy in re-writing a memory. e A sample is first positively trained by σ+ applied for 250 s. We then estimate the time dependent power dissipation (P+−) when a stress σ− (|σ+| = |σ−|) is applied in the negative direction on this sample. We normalize P+− by the power dissipation (P++) when σ+ is applied on the initially positively trained sample. Each curve (solid circles) is obtained for one pair of σ+ and σ− values. The 3-D plot is generated by varying the magnitude of such pair [σ+, σ−].
stiffness versus strain curve highly asymmetrical.
This collective nonlinearity, i.e. nonlinearity emerging from the interaction of many filament, is the dominant effect for small and moderate strains 22 (Fig.4a). At larger strains, the constitutive nonlinear stiffening of individual cross-linkers contributes as well 21 (see Fig.S4). This leads to a strong strain stiffening of the network at larger values of γ, where the collective nonlinearity combines with the constitutive nonlinearity to produce additional asymmetric strain stiffening.
In order to compare our model to the experiment, we estimate the degree of nematic order produced by training. Here, we do not consider in the simulation the training and cross-linking dynamics. Rather, we rely on the observed post-training residual strain γ r (Fig.3b, also see Fig.S5), which provides a direct measure of the change in the elastic reference state. Presuming an affine deformation, we then calculate the expected nematic order induced by filament reorientation by this shear (see Method). Using this simple, geometric model to infer the degree of nematic order S(γ r ), we numerically construct an equivalent semiflexible network for each The nematic order encodes the history dependent shear response of actin networks. a, Measurement and prediction of the relative stiffening as a function of shear strain γ, for different training times T+ as described in the text. The training stress σ+ = 4 Pa in all cases. Full markers represent experimental data while continuous curves the corresponding predicted data from simulations. b, Cartoon depicting that for a positively trained sample, positive RO strains put most filaments are under tension which lead to a larger overall network stiffness. Conversely, negative strains put most filaments under compression, the network enters a soft bending mode and the overall stiffness decreases. c, The residual strain measured experimentally after training is used to compute an equivalent nematic order S under the assumption of affine deformation. This nonlinear behavior is magnified as the nematic order increases. d, We observe an excellent agreement between the relative stiffening measured experimentally (discrete points) at |γ|=0.4 for different values of residual strains (e.g. different training time) and predicted by the simulations (red line with symbol) for the corresponding nematic order, as well the prediction from the simple mean field model (black dashed line).
The tangent modulus K(γ), normalized with respect to its value before training is shown for one untrained (black) network and two networks trained for 200s (blue) and 400s (red) as lines in Fig.4a. We find excellent agreement with the experimental data. We then plot ∆K as a function of residual strain (Fig.4d). The agreement between the observed and simulated mechanomemory development strongly suggests that the development of nematic order underlies the observed changes in elastic nonlinearities. Interestingly, we see evidence of such nematic order in the network under shear deformations (Fig.S6), which is qualitatively consistent with our model assumptions.
To further investigate the role of nematic order in the network and consequently the emergence of buckling in its nonlinear response, we build a mean field network model composed of linear elastic filaments with frozen nematic order S, and incompliant cross-linkers. We account for filament buckling simply by neglecting the mechanical contribution from filaments under compressive loads exceeding their Euler criterion 24 (see Supplement). Figure 4d plots ∆K predicted by the mean-field model (dashed black line) agrees with both the simulation and experiment. Based on the mean field model, we suggest that proposed dependence of filament buckling on the degree nematic alignment in the gel serves as a minimal model of mechanomemory effects.
Our study demonstrates the reversible formation of mechanomemories in physiologically relevant cross-linked actin networks. Previous studies of strain-history dependence in the rheology of biopolymer networks observed irreversible effects attributed to permanent changes in network bundling 7, 14,15 . Importantly, using the identical cyclic hardening protocol used in 7 , we do not observe work hardening (Fig.S7). This confirms the effect we describe here is distinct from irreversible work hardening observed in previous studies. Future work is needed to understand how training stress and duration as well as network composition can control the formation of either reversible or irreversible stress-dependent changes in the mechanical response.
The field-dependent hysteresis in mechanical response we describe is highly reminiscent of other hysteresis phenomena (e.g. magnetism). For instance, our data suggests that the stress-dependent changes in mechanical response arise from the formation of nonequilibrium, but long-lived domains of partial nematic order in the network that strongly modifies the nonlinear response of the network to strain. Care must be taken to distinguish such hysteresis from the ubiquitous rate-dependent hysteresis observed in viscoelastic materials. Our study demonstrates that that the stress-dependent changes in mechanical response can be understood solely in terms of the development of long-lived nematic order.
Previous work on nematic elastomers 25,26 and wormlike micelles 27 has demonstrated the importance of nematic order on mechanical response. Here, the ability to rewrite and erase memories comes from the fact that the lifetime of the nematic order is highly stress dependent. Networks of elastically nonlinear elements (e.g., filaments subject to softening under compression in response to buckling) that can structurally reorganize under load (develop nematic order) point to an avenue for the development of smart adaptive materials based on these cytoskeletal structures. There remain open questions regarding how mechanomemory depends on network density and composition. We have focused here on actin networks cross-linked with the flexible and dynamic cross-linker filamin. We observe qualitatively similar, but quantitatively distinct, behavior using the cross-linking protein α-Actinin which is more rigid and shorter, but with similar binding kinetics (Fig.S8). This suggests that cross-linker properties may prove critical to controlling the details of the stress-dependent changes in nematic order. Moreover, understanding the role of network density and, specifically on the transition from non-affine to affine deformations 22,28,29 is also likely to be important. Exploring these properties in biopolymer networks will enable a more complete understanding of the control parameters to inform design of materials constructed from synthetic analogs.
Our study demonstrates the possibility to reversibly alter the mechanical response of network to enable one to write, read and erase mechanomemories. In fact, this system may have the potential to encode multiple memories 30 through the formation of multiple nematic domains with differing nematic directors, leading perhaps to new classes materials with a programmable and complex nonlinear elastic response. Finally, our results also suggest a mechanism of observed force adaptation in living cells and tissue 31 .

Protein preparation
Monomeric actin (G-actin) was purified using protocol adapted from Ref. 32 from rabbit skeletal muscle acetone powder (Pel Freeze Biologicals, Product code: 41008-3). Small aliquots of actin solution were drop frozen in liquid nitrogen, and stored at -80C. For each experiment fresh aliquots of frozen actin were thawed and used. FLN was purified from chicken gizzard and aliquots were also stored at -80C.

Bulk rheology
All rheological measurements were performed on a Bohlin Gemini HR Nano (Malvern Instruments) rheometer at a temperature of 25 0 C, using a 40 mm diameter acrylic plate and 160 µm gap. The acrylic plate has much lower moment of inertia compared to a similar metal plate and hence rheological measurements are more robust under fast direction switching of the plate rotation. We confirmed results were independent of plate geometry (Fig.S9). We also used a humidity chamber sealed with vacuum grease to reduce water evaporation from the sample during the rheology experiments and observed that the sample's rheological properties were essentially unchanged for periods up to ∼ 4 hours after loading the sample in the rheometer. The mechanomemory effect described in this manuscript is a robust effect, observed 100% of the time in > 35 independent samples.
After loading the sample, the build-up of linear elastic G and viscous G moduli was measured as a function of time at a fixed frequency f = 0.5 Hz with a very small strain amplitude γ 0 = 0.02 (ensuring linear response) to monitor the polymerization process indicated by the increase in magnitude of G and G with time. After approximately 1 hour, the moduli saturate indicating a fully polymerized starting state of cross-linked actin network. This state is the untrained state of the network. The tangent shear modulus K(γ) is measured from the Lissajous plots (Fig.1d) by numerically calculating δσ δγ | γ along the increasing strain magnitude under various training conditions. Memory is quantified by the asymmetry of tangent shear modulus ∆K in positive and negative direction at strain |γ| = 0.4, i.e. ∆K(|γ|) = K(+γ) -K(−γ). The choice of this strain controls the absolute value of the strength of memory, but does not affect changes in memory due to training or our analysis. The biggest source of error in shear modulus K (or ∆K) comes from the noise due to differentiation of high resolution data (we get K by differentiating the stress with respect to strain). Here, we first numerically fit a spline smoothing curve (averaged over 10 adjacent points). Then we take the difference of this smooth curve from the raw data. The standard deviation (SD) of such differences gives the error in K. The error bars in ∆K are twice the standard deviations, since the uncertainty gets doubled when taking the difference.
To get the phase information, the complex tangent shear modulus K * (σ, ω) was measured by superposing a small sinusoidal component of amplitude δσ and angular frequency ω on a constant background stress σ with δσ = 0.1σ. Then we measured the amplitude of the sinusoidal strain δγ at the same angular frequency ω and the phase difference ∆ between the sinusoidal stress and the strain. The magnitude of the tangent shear modulus is given by K(σ, ω) = δσ δγ | σ . Under steady shear, P = σγ gives the instantaneous injected power per unit volume. However, since the system is visco-elastic, some part of this injected energy is dissipated and some part is stored. When the training stress is switched off, the system can relax back partially using the stored energy. We can also get insight about the dissipation in the system by superposing a small amplitude sinusoidal component on top of the large D.C. stress component. By this method we can measure the phase difference between the applied sinusoidal stress and resulting sinusoidal strain. When we plot the sinusoidal strain vs sinusoidal stress over a complete cycle, we get a closed curve (Lissajous plot) that can be well approximated by an ellipse in our case. The area under the curve is proportional to the total dissipated energy over a complete cycle. Here, it should be clear that the dissipated energy will also depend on the amplitude and frequency of the applied sinusoidal signals. Thus, by this method we cannot estimate the absolute power dissipation. How-ever, if we use the same amplitude and frequency of the input sinusoidal component, we can compare how the power dissipation occurs as a function of time during the creep behavior under D.C. stresses. The energy dissipated over a complete cycle is E cycl = π δσ δγ sin∆, leading to an average instantaneous power dissipation of P (t) = (ω/2π)E cycl . We define P αβ (t) to be the average instantaneous power dissipation in the β direction in a sample trained in the α direction. Here α and β can assume only the values ±. Under steady shear, we calculate the injected power and subtract out the recoverable power at every instant to estimate the time dependent power dissipation (Fig.S2). However, to do this one has to assume that the injected power at the start up of creep response corresponds to the recoverable power at the start of recovery curve.

Confocal microscopy
20% molar ratio of labeled actin monomers (actin from rabbit skeletal muscle conjugated with Alexa fluorophore 568, Life Technologies Corp., product part no. A12374) was added to 80% G-actin and then the combination is mixed gently with 1X Ca-G-buffer. Ca-G-buffer containing both labeled and unlabeled G-actin and 10X actin polymerization buffer containing FLN were mixed gently for 10 s just before loading the sample into a home built glass chamber (thickness ∼ 100µm) mounted on the sample stage of a laser scanning confocal micro-scope (Nikon) with a 60X water immersion objective. Z-stack images were recorded approximately after 1 hour of polymerization.

MODEL
We developed a mathematical model of the system as a two-dimensional network of N semiflexible filaments with identical elastic properties and length L. To construct the numerical network, each filament is placed at random in a box of area A with an orientation picked from a distribution P (θ), where θ denotes the angle between the filament and the nematic director n of the network.

Single filament and crosslink mechanics
Using classical Bernoulli-Euler beam theory, the mechanical properties of the filaments are characterized by the linear stretching and bending moduli, µ and κ respectively. The elastic length scale or bending length λ b = κ/µ determines the relative contribution of the filament's bending and stretching modes to the total energy. The thermal persistence length of the filaments is given by = κ/k B T . The energy of a single filament can then be written as where we used an arclength parametrization of the filament with r(s) the position of the filament cross-section after deformation, and θ the rotation of the filament caused by that deformation.
We implement a mechanical model for the filamin cross links that is based on previous experimental measurements 33 . The filamin contour length is ∼100 nm, but is comprised of subunits of rest length ∼10 nm that unfold under force. Previous experimental results show that, during unfolding, the subunits follow a worm-like chain model 33 . In the simulation, the saw-tooth forceextension curve is reduced to two regimes, reflecting the regime of subsequent unfolding and the nonlinear stiffening observed during stretching of the last segment. This was adopted in order to reduce problems arising from numerical instability. For lengths below 150 nm, the 1d spring constant is µ lc = 100pN . In the stiffening regime above 150nm, there is a linear increase in the spring constant of 10 pN/nm for further extension.
The geometry of the network is characterized by the scalar amplitude of the nematic order parameter S, found from the width of the distribution P (θ): (2) Our results are remarkably insensitive to the exact form of distribution of filament angles (Fig.S10).
We have been using the following distribution in our simulations: where I 0 (α) the modified Bessel function of the first king of order 0. The parameter α is tied to the nematic angle and is found numerically by solving S = I 1 (α)/I 0 (α).
Another geometric characteristic of the network is the mean distance between consecutive crosslinks along a given filament, l c , which is inversely proportional to the network density and depends on the nematic order parameter S. This quantity is necessary to assess the propensity for buckling. Long spans between consecutive cross links buckle at a lower compressive stress than shorter ones. This quantity may be computed by considering a filament lying at an angle θ with respect to the nematic director. The probability of a second filament crossing at an angle ψ with respect to the first is given by L 2 |sin ψ| /A. In order to obtain the crossing probability independent of the relative angle ψ, we need to integrate the above quantity over ψ with a weight P (θ + ψ). Doing this gives: This quantity is simply the probability that a given filament oriented with angle θ to the nematic direction is crossed by a second filament with its direction chosen from the angular probability distribution and its position chosen randomly. We will again assume that our system contains a large number of filaments N in a large area A, so that P cross (θ) is small; in this limit, the probability distribution p n (θ) for the number of crosslinks on a filament oriented with angle θ to the nematic direction can be approximated by the following Poisson distribution: where n c (θ) is the mean number of cross links on a filament oriented with angle θ, and is written: with ρ = N/A the filament density. We then calculate the distance Λ c between the two most distance crosslinks on a filament and the average number of intervals between two crosslinks on a filament N c . These quantities are given by: Finally, the ratio of Λ c and N c yields the mean distance between crosslinks: where · denotes the averaging operation over θ with weight P (θ).

Simulation details
We consider a system composed of N semiflexible filaments connected together by cross links at points of intersection. The cross links are extensible nonlinear elastic elements as described above, but they supply no constraint torques, being free to rotate. The filaments are spatially discretized by placing nodes at each cross link, as well as at regular intervals between adjacent cross links, which ensures a sufficient level of spatial refinement for numerical convergence. Typically, networks are generated with 4000 filaments discretized into 100 000 finite elements. The construction of each network is based ... on the two geometric parameters defined above: (a) the nematic order parameter S and (b) the mean filament density l c . To build the network, we add filaments sequentially to the simulation box with typical dimensions of A = 64L 2 . We numerically evaluate the strain energy using the discretized version of equation (1): where d and d 0 are respectively the current and original length of the segments, and l 0 is the average rest length of two adjacent segments forming an angle β. The nematic order parameter is set by choosing the direction of the filaments from a Gaussian distribution with a width determined to give the desired nematic order parameter S. The centers of the mass of the filaments remain uniformly distributed during the construction of all networks. We have explored the mechanics of networks constructed using other anisotropic orientational distributions of the filaments, and found only small changes in the collective mechanics of the networks when both the network density and nematic order parameter were held constant. See Fig.S10.
We use the Lees-Edward method 34 to simulate the shearing deformation of the network with periodic boundary conditions, and the static equilibrium of the system is obtained by energy relaxation using a quasi-Newton minimization algorithm 35,36 .

From residual strain to nematic order
The residual strain measured after training implies a realignment of the filaments in the network. Assuming that the residual strain results from the affine deformation of a formerly isotropic network, we compute the nematic order parameter S induced by training. Choosing an arbitrary filament from an isotropic 2D network of filaments (length L), we write its initial orientation as with θ the angle made by the filament with respect to the x axis. We apply a simple shear along direction x in the xy plane. Under this deformation the filament rotates make the angle θ with respect to the same x axis, given by The new angle θ is thus θ = arctan sin θ cos θ + γ sin θ .
Considering the ensemble of filament orientations given by the unit vectors u = (cos θ, sin θ) in the isotropic network, we apply the mapping Eq. 13 to obtain the new (anisotropic) distribution u = (cos θ , sin θ ) of filament orientations after the shear deformation λ = δ + γxy. From this distribution we directly compute the resulting nematic order parameter tensor by averaging over the initially isotropic distirbution of angles θ. Denoting this average by . , we write the order parameter S ij = u i u j − 1 2 δ ij as: Finally, after diagonalization, we find the magnitude of the order parameter S to be given by This provides an explicit expression for the nematic order as a function of the residual strain γ.

Mean field approach
We develop a mean field calculation of the shear modulus as a function of nematic order. We shear the network along the +x direction, where the nematic director makes an angle φ with respect to x (Fig.M1b). We assume that filaments buckle only when under a compressive load larger than the Euler buckling load p c , which depends on the length of the filament segment between consecutive cross links and on the angle that that filament makes with the shearing direction. We further assume that one may neglect the (small) amount of elastic energy stored in filaments post-buckling. Under these assumptions, a filament oriented making an angle ψ with the x-axis has an energy per unit length if not buckled. Based on the Euler buckling criterion, there is a wedge of width 2ω Euler ( ) around the compression direction of the shear so that filament segments of length greater than oriented within this wedge will buckle. We wish to exclude such buckled filaments from the computation of the elastic energy stored in the network using Eq. 16. We replace this segment-length dependent quantity 2ω Euler ( ) with an averaged excluded wedge angle, θ c . This is determined self-consistently by setting the average filament segment length (averaged only over the wedge of half angle θ c to be equal to the load-dependent critical length for Euler buckling l buckle (γ). This results in a self-consistent integral equation to be solved for θ c , and which depends on both the degree of nematic order and shear strain γ: where · notes denotes the average over the wedge for a given degree of nematic order P (θ). This is solved numerically. We then compute the total network energy by integrating Eq. 16 over the filament orientation and length 23 , excluding wedge of half angle θ c surrounding the compression direction, oriented along θ comp = -45 degrees w.r.t. the shearing direction. The total energy of the network is Using the fact that E(φ) = K(φ)γ 2 /2 for small shear, we write the tangent modulus as θcomp−θc 0 dψ cos 2 ψ sin 2 ψ/2 × ∞ n=2 L(n − 1)p n (ψ − φ) n + 1 + 2µρ π θcomp+θc dψ cos 2 ψ sin 2 ψ/2 The mean-field value of ∆K computed using Eq. 19 is shown in Fig. 4d as a dashed black line. The mean-field model is consistent with both the simulation (red line and symbols) and the experiment (points).   Fig.2d in the main text. This plot indicates that, starting with a highly asymmetric state (large memory), it is possible to bring back the sample to a state that is close to the untrained (no memory) state by applying suitable opposing stress pulses. After waiting for 300 s, a stress of magnitude 5 Pa is again is applied for 300 s and then the stress is switched off (σ = 0 Pa) for 300 s. We estimate the shear rate as a function of time by numerically differentiating the strain data during both creep (γc(t)) and recovery (γr(t)) as indicated in Fig.3 in the main text. The difference shear rate is measured by ∆γ(t) =γc(t) −γr(t). The quantity σ ∆γ(t), gives the instantaneous power dissipation. P++(t) and P+−(t) represent the power dissipation when the creep stress direction is same or opposite to the training direction, respectively. b, P++(t) and P+−(t) as a function of time. The integrated area between P+−(t) and P++(t) indicated by the gray region gives the total dissipated energy ∆ EM in re-writing a memory. The symbols represent the raw data and the lines represent the spline curve (obtained by 5 points moving averaging). c, ∆ EM as a function of D.C. stress pulse magnitude. The error bars represent the standard deviation of the fluctuations in power vs the time data in (c).   represents the corresponding FFT image where the colors from blue to red to bright yellow represent increasing magnitude of Fourier amplitudes. We see a frozen in nematic order in the direction of applied shear (shown by an arrow), resulting from the residual strain. FIG. S10. Linear modulus versus strain for nematic networks prepared with different forms of their orientational order. The (red) uniform networks have step orientational distributions so that all angles within ∆θ = π/3 the nematic director are equally probable and all other orientations are disallowed. The (blue) Gaussian distributions (σ 2 = π/8) replace the step distribution with a Gaussian one, as used in all simulations reported in the manuscript. The comparisons are done for two different samples (one pair of curves for each) to also study the sample to sample variations. The difference in modulus between these networks is small in all cases.