A minimal rupture cascade model for living cell plasticity

CML arises from a hematopoietic stem cell transformation following the formation of the BCR-ABL oncogene by a single reciprocal chromosomal translocation t(9;22) [40]. In CML, BCR-ABL⁺ cells of the myeloid lineage proliferate uncontrollably, the bone marrow density increases considerably, and their mechanical properties change during disease progression [41]. In transformed cells [42], BCR-ABL was shown to bind actin filaments, to inhibit their polymerization and to disorganize the CSK into punctuate, juxtanuclear aggregates [43–45]. We used AFM to indent single hematopoietic purified CD34⁺ cells from CP-CML patients at diagnosis and healthy donor bone marrows [36] (see appendix A). We indented the cells by moving vertically the AFM cantilever (figure 1(a)) toward their perinuclear central region (figure 1(b)), at constant speed V₀ = 1 μm s⁻¹ until a cantilever set-point force was reached. Then the cantilever was withdrawn from the sample at the same constant speed (−V₀) back to its starting position (figures 1(c) and (d)). From these indentation experiments, the cell shear relaxation modulus G(Z) was estimated as the second-order derivative of the FIC (see equation (B13) in appendix B) [36, 39]. From the approach and retract FICs, the initial G_i and global G_g shear moduli, and the dissipation loss D_l were retrieved (see figure S1 is available online at stacks.iop.org/NJP/20/053057/mmedia and supplemental material¹¹ ). When these FICs do not superimpose (figures 1(c) and (d)), D_l quantifies the percentage of work not restituted during retract and dissipated partly as viscous loss [31, 36]. Interestingly, the experimental FICs display fluctuations that locally exceed the background thermal fluctuations of the AFM cantilever [36] (figures 1(c) and (d)). We developed a wavelet-based detection method of these singular events that amounts to detect local curvature minima G_m and maxima G_M in the FICs (figures 1(e) and (f)) (see also appendix B, and figures S2 and S3 (see footnote 11)). For each disruption event, we computed the force drop ΔF, the penetration length ΔZ, and the released energy E = ΔF.ΔZ (figure 1(e)). FICs without disruption events were not included in the statistics.

Zoom In Zoom Out Reset image size

We collected two large sets of single cell FICs from 5 CP-CML patients (1301 FICs—49 cells) and 5 healthy donors (1671 FICs—60 cells). We detected 6161 singular rupture events distributed on 1153 FICs in CP-CML cells as compared to 6765 rupture events distributed on 1111 FICs in normal cells. Thus only 11.4% (148/1301) of FICs do not display rupture events for CML cells which is significantly lower than 33.5% (560/1671) for normal cells. The computation of the normalized histograms (p.d.f. for probability density function) of ${\rm{\Delta }}Z$, ΔF, and E of these rupture events, revealed rather wide distributions with fat-tail (see figure S4 (see footnote 11)). When focusing on ΔZ (nm), which can also be interpreted as the time duration Δt = ΔZ/V₀ (ms) of the rupture event, we got very satisfactory fits of the p.d.f. of ${\mathrm{log}}_{10}({\rm{\Delta }}Z)$ with the sum of two distinct Gaussian distributions (see equation (B14) in appendix B) for both normal and CML cells (figures 2(a) and (b)). Thus, we identified two populations of rupture events, a subpopulation 1 of rather short duration ($\overline{{\rm{\Delta }}{t}_{1}}\sim 30\,\,\mathrm{ms}$) and weakly penetrating ($\overline{{\rm{\Delta }}{Z}_{1}}\sim \,30\,\mathrm{nm}$) failures, and a subpopulation 2 of longer ($\overline{{\rm{\Delta }}{t}_{2}}\sim \,50\,\,\mathrm{ms}$) and deeply penetrating ($\overline{{\rm{\Delta }}{Z}_{2}}\sim \,50\,\mathrm{nm}$) failures. But what distinguishes CML from normal cells is the higher percentage (α = 0.34 compared to 0.26) of the larger rupture events (subpopulation 2) in the CP-CML cells, as an indication of their greater mechanical brittleness.

Zoom In Zoom Out Reset image size

When investigating the correlations between ΔZ, ΔF and E, we found a rather strong correlation between ${\mathrm{log}}_{10}({\rm{\Delta }}F)$ and ${\mathrm{log}}_{10}({\rm{\Delta }}Z)$ for both the CML (Pearson's correlation coefficient r = 0.68) and the normal (r = 0.62) cells. Log₁₀ (E) and ${\mathrm{log}}_{10}({\rm{\Delta }}Z)$ were also found significantly correlated in CML (r = 0.83) and in normal (r = 0.80) cells. These correlations issue from the existence of two clouds of points in the respective scatter plots (figures 2(c) and (d)), corresponding to subpopulation 1 of small indentation depth, short duration, small force drop and low released energy failures, and to subpopulation 2 of large indentation depth, long duration, large force drop and high released energy failures. We will classify the former as ductile and the latter as brittle rupture events.

The pertinence of the fitting of the p.d.f. of ${\mathrm{log}}_{10}(E)$ (figure 2(e)) by the sum of 2 Gaussians is compelling when using a logarithmic representation (figure 2(f)). For the CML cells, when fixing the relative percentages of ductile (1 − α = 0.66) and brittle (α = 0.34) rupture events as previously estimated, we obtained the parameter values reported in table 1 with the following arithmetic and geometric means $\overline{{E}_{1}}=216$ k_BT (resp. $\overline{{E}_{2}}=1547$ k_BT), and ${\tilde{E}}_{1}=141$ k_BT (resp. ${\tilde{E}}_{2}=776$ k_BT). For normal cells, when fixing α = 0.26 (1 − α = 0.74), we ended with consistent parameter values (table 1) with $\overline{{E}_{1}}=188$ k_BT (resp. $\overline{{E}_{2}}=1158$ k_BT), and ${\tilde{E}}_{1}=138$ k_BT (resp. ${\tilde{E}}_{2}=741$ k_BT). Note that the standard deviations of the Gaussian distributions for ${\mathrm{log}}_{10}(E)$ are slightly larger for CML than for normal cells and this for both ductile and brittle events.

Table 1.  Distributions of energy released during ductile and brittle rupture events in normal and CML cells: simulations versus experiments. Characteristics ($\overline{N}$, σ_N, $\mu =\overline{{\mathrm{log}}_{10}(E)}$, $\sigma ={\sigma }_{{\mathrm{log}}_{10}E}$, $\overline{E}$, $\tilde{E}$) of the numerical cascades simulated with equation (2) for parameters values a₀, $\hat{a}$, ΔE₀ = 12 k_BT, ${\rm{\Delta }}{E}^{* }=4$ k_BT versus experimental data (figures 2(e) and 2(f)).

	Normal Cells	CML Cells
Ductile
Model parameters	a0 = 1.3, $\hat{a}=15$, Δa = 0.39	a0 = 1.3, $\hat{a}=15$, Δa = 0.48
Simulations	$\overline{N}=9$, σN = 3	$\overline{N}=8$, σN = 4
	μ1 = 2.17, σ1 = 0.31	μ1 = 2.13, σ1 = 0.38
	${\overline{E}}_{1}=191\,$ kBT, ${\tilde{E}}_{1}=148\,$ kBT	${\overline{E}}_{1}=198\,$ kBT, ${\tilde{E}}_{1}=135\,$ kBT
Experiments	μ1 = 2.14, σ1 = 0.34	μ1 = 2.15, σ1 = 0.40
	${\overline{E}}_{1}=188\,$ kBT, ${\tilde{E}}_{1}=138\,$ kBT	${\overline{E}}_{1}=216\,$ kBT, ${\tilde{E}}_{1}=141\,$ kBT
Brittle
Model parameters	a0 = 1.3, $\hat{a}=46$, Δa = 0.35	a0 = 1.3, $\hat{a}=53$, Δa = 0.38
Simulations	$\overline{N}=21$, σN = 7	$\overline{N}=22$, σN = 8
	μ2 = 2.84, σ2 = 0.42	μ2 = 2.92, σ2 = 0.47
	${\overline{E}}_{2}=1104\,$ kBT, ${\tilde{E}}_{2}=692\,$ kBT	${\overline{E}}_{2}=1494\,$ kBT, ${\tilde{E}}_{2}=832\,$ kBT
Experiments	μ2 = 2.87, σ2 = 0.41	μ2 = 2.89, σ2 = 0.51
	${\overline{E}}_{2}=1158\,$ kBT, ${\tilde{E}}_{2}=741\,$ kBT	${\overline{E}}_{2}=1547$ kBT, ${\tilde{E}}_{2}=776\,$ kBT

When comparing the mean released energies during ductile (∼200 k_BT) and brittle (∼1300 k_BT) rupture events to the dissociation energies of α-actinin/actin binding (∼4.3 k_BT) and of filamin/actin binding (∼3.6 k_BT) [46], we obtain rough estimates for the number of ABP unbindings during ductile (∼50) and brittle (∼325) rupture events. Thus, more than 6 times energy is released during brittle failures that last (∼50 ms) not more than twice the duration (∼30 ms) of ductile failures, meaning that the mean rate of released energy is significantly higher during brittle failures. Interestingly, the computation of the global G_g and initial G_i shear moduli (see figure S1 (see footnote 11)) confirms that CML cells are stiffer ($\overline{{G}_{g}}=1.13\pm 0.31$ kPa, $\overline{{G}_{i}}=0.77\pm 0.26$ kPa) than normal ones ($\overline{{G}_{g}}=0.42\pm 0.05$ kPa, $\overline{{G}_{i}}=0.28\pm 0.05$ kPa) [36]. In particular, about 37% of FICs bear mean shear relaxation moduli (G_g) larger than 1 kPa in CML cells with very low (7%) counter part in normal cells (see figure S5 (see footnote 11)). Altogether, these observations show that CML cells display a significant proportion of highly tensed perinuclear zones propitious to localized brittle failures by disruption of cross-linked CSK domains impeding complete shape recovery after deformation. This interpretation is strengthened by the experimental observation of an important structural alteration of the actin CSK of immature TF1 cells consecutive to BCR-ABL oncogene transduction. In particular, confocal fluorescence microscopy revealed that in this model of human CML [40–45], juxtanuclear actin aggregates were found in almost 30% of the BCR-ABL-transduced TF1 cells at the expense of the cortical F-actin [43, 45] (figure 3(b)). Very likely, these solid structures were induced by the oncogene since they were rarely observed in the parental TF1 cell line (figure 3(a)).

Zoom In Zoom Out Reset image size

P.d.f.s with heavy tails have been widely observed in various domains of fundamental and applied sciences [50]. The most popular fat-tail distributions are power-law and log-normal distributions that often have been considered as competing models of experimental data [48, 49]. Power-law distributions are commonly thought to be a statistical characteristics of systems that display space and/or time scale invariance properties [50] with as notable examples scale-free networks [51] and self-organized critical systems [52–54]. Log-normal distributions are paradigmatic fat-tail distributions generated by self-similar fragmentation and more generally multiplicative processes [55, 56] with as historical example the energy cascade model of fully developed turbulence [57]. But, as originally proposed by Kesten [47], power-law and log-normal distributions are indeed intrinsically connected when combining multiplicative and additive random processes:

$$\begin{eqnarray}&&{S}_{t}={a}_{t}{S}_{t-1}+{b}_{t},\end{eqnarray} \tag{ 1 }$$

where S_t is the size of a population at time t, a_t the random positive growth factor, and b_t a small positive random increment. In the absence of the additive b_t term, one recovers the multiplicative Gibrat's law [58] of proportional growth which is nothing but a random walk in log-size leading to a log-normal distribution of S_t. But this distribution is not stable since when the time increases the process S_t either asymptotically shrinks stochastically to zero ($\overline{\mathrm{ln}({a}_{t})}\lt 0$) or diverges to infinity ($\overline{\mathrm{ln}({a}_{t}}\gt 0$). Kesten [47] showed that provided ($\overline{\mathrm{ln}({a}_{t})}\lt 0$), adding the random variable b_t prevents the divergence of the log-normal distribution of S_t which progressively switches and converges to a stationary power-law distribution with exponent α given by the strictly positive solution of the equation $\overline{{a}_{t}^{\alpha }}=1$.

Our model of CSK cascading rupture events is deliberately reductionist with a minimal number of parameters. At the cascade step n + 1, the energy released ΔE_n+1 satisfies a Gibrat's multiplicative law [58], meaning that it can be expressed as a percentage of the energy ΔE_n released at the previous cascade step:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{n+1}={a}_{n}{\rm{\Delta }}{E}_{n},\,\,\,{\rm{with}}\,\,{a}_{n}\sim { \mathcal N }({a}_{0}{{\rm{e}}}^{-n/\hat{a}},{\rm{\Delta }}{a}^{2}),\end{eqnarray} \tag{ 2 }$$

where a_n is a Gaussian random variable with initial value a₀ > 1, whose mean $\overline{{a}_{n}}$ decreases exponentially (CSK structural relaxation) with a characteristic depth $\hat{a}$ and a fixed standard deviation ${\sigma }_{{a}_{n}}={\rm{\Delta }}a$. Note that the relationship between the cascade index n and time is not defined in our model since there is no reason a priori to assume that the CSK rupture cascade steps occur at regular time intervals. Starting from some initiation threshold ΔE₀, ΔE_n will asymptotically tend to zero because of the exponentially decreasing multiplicative constant a_n which will become smaller than 1 after some finite number of steps. The cascade will end when ${\rm{\Delta }}{E}_{N+1}\lt {\rm{\Delta }}{E}^{* }$, where ${\rm{\Delta }}{E}^{* }\geqslant 0$ is a cascade arrest energy cut-off. During the N steps of the CSK rupture cascade, the total released energy is:

$$\begin{eqnarray}&&E=\displaystyle \sum _{n=1}^{N}{\rm{\Delta }}{E}_{n}.\end{eqnarray} \tag{ 3 }$$

Our random cascade model mainly depends on three parameters: a₀, $\hat{a}$ and Δa. N is a random variable that depends on the realization of the cascade. Its p.d.f. P(N) is well approximated by a Gaussian distribution (figure 4). To account for the experimental data, all our simulations of the released energy cascade model defined by equation (2) were performed for rather small values of Δa ≪ 1. In this limit, E (equation (3)) can be approximated by:

$$\begin{eqnarray}&&E\simeq {E}_{0}\displaystyle \sum _{n=1}^{N}{a}_{0}^{n}\,{{\rm{\Pi }}}_{i=1}^{n}{{\rm{e}}}^{-i/\hat{a}}\,={E}_{0}\displaystyle \sum _{n=1}^{N}{a}_{0}^{n}\,{{\rm{e}}}^{-{\sum }_{i=1}^{n}i/\hat{a}}.\end{eqnarray} \tag{ 4 }$$

From the well known identity ${\sum }_{i=1}^{n}i=n(n+1)/2$, we get

$$\begin{eqnarray}&&E\simeq {E}_{0}\displaystyle \sum _{n=1}^{N}{a}_{0}^{n}\,{{\rm{e}}}^{-\tfrac{n(n+1)}{2\hat{a}}}.\end{eqnarray} \tag{ 5 }$$

Thus a good approximation of the mean of E can be obtained by replacing N by the nearest integer $[\overline{N}]$ of its mean $\overline{N}$ in equation (5):

$$\begin{eqnarray}&&\overline{E}\simeq {E}_{0}\displaystyle \sum _{n=1}^{[\overline{N}]}{a}_{0}^{n}\,{{\rm{e}}}^{-\tfrac{n(n+1)}{2\hat{a}}}.\end{eqnarray} \tag{ 6 }$$

We recall that the p.d.f. of the random variable N was shown to be well approximated by a Gaussian (figure 4). Now from the Taylor expansion of $\mathrm{ln}(E)=\mathrm{ln}(\overline{E})+(S-\overline{S})/\overline{S}$ at first order, we can show by taking the ensemble average on both sides of this equation that $\mathrm{ln}(\overline{E})$ is well approximated by

$$\begin{eqnarray}&&\overline{\mathrm{ln}(E)}\approx \mathrm{ln}(\overline{E}).\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

We have confirmed numerically the pertinence of this approximation in all our simulations with a rather good accuracy (<5% error). Now, if we assume that the energy cascade ends when the difference between $\overline{{\rm{\Delta }}{E}_{n}}$ and the threshold ΔE* is of the order of Δa, i.e. ${{a}{\rm{e}}}^{-\overline{N}/\hat{a}}\simeq K{\rm{\Delta }}a$, then

$$\begin{eqnarray}&&\overline{N}\simeq \hat{a}\,\mathrm{ln}\left(\displaystyle \frac{{a}_{0}}{K{\rm{\Delta }}a}\right),\end{eqnarray} \tag{ 8 }$$

where K is a constant that has been numerically estimated of order 1 (e.g. K = 1.52 ± 0.05 for a₀ = 1.3, $\hat{a}=45$ and Δa = 0.36, see figure 5). Equation (8) shows that the mean size $\overline{N}$ of the energy cascade increases as expected when increasing a₀ and $\hat{a}$, but decreases when increasing Δa, since the probability to have a ${\rm{\Delta }}{E}_{n}\lt {\rm{\Delta }}{E}^{* }$ after a limited number n of cascade steps increases (see figure S6 (see footnote 11)).

Zoom In Zoom Out Reset image size

In our simulations, we fixed a₀ = 1.3 and started iterating equation (2) with the initial threshold ΔE₀ = 12 k_BT. Note that ΔE₀ is a multiplicative factor that only shifts the p.d.f. of ${\mathrm{log}}_{10}(E)$ without affecting its shape. When using the cut-off ${\rm{\Delta }}{E}^{* }=0$ (${\rm{\Delta }}{E}_{n}\gt 0,\forall n$ ), after some transient (n ≳ 5), the p.d.f. of ΔE_n obtained from 1.2 × 10⁶ realizations with parameter values $\hat{a}=45$ and Δa = 0.36 is well approximated by a log-normal Gibrat's law [58], as expected from multiplicative process (figures 5(a) and (b)). What is more surprising is the fact that the sum E of these log-normal variables turns out also to be well approximated by a log-normal distribution and not by a normal distribution as expected from the central limit theorem. Indeed, this theorem does not apply since the random variables ΔE_n are not only not identically distributed but correlated. Beaulieu [59] has recently proved an extended log-normal limit theorem that states that the limit distribution of the sum of nonidentically distributed correlated log-normal random variables is a log-normal distribution. When investigating how the energy ${E}_{{n}^{* }}={\sum }_{n=1}^{{n}^{* }}{\rm{\Delta }}{E}_{n}$, released during the first n^* steps of the cascade converges to the total released energy E_N, we confirmed a rather fast convergence to an asymptotic log-normal distribution (figures 5(c) and (d)). When numerically simulating surrogated uncorrelated released energy cascades where at each step, ΔE_n is randomly generated with the p.d.f. P(ΔE_n) (figure 5(a)) but independently of the previous cascade steps, then the sum ${E}_{{n}^{* }}$ no longer converges to a log-normal distribution but to a normal distribution as expected for uncorrelated random variables (figures 5(e) and (f)). This theoretical argument suggests that the log-normal distribution observed experimentally (figures 2(e) and (f)) are evidence that the different cascade steps of ABP unbinding are correlated. This theoretical understanding applies to all the numerical simulations reported in this study, for ductile and brittle rupture events, in normal as well as in CML cells (table 1, figures 6 and 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We performed numerical simulations of the released energy cascade model to reproduce quantitatively the p.d.fs of E obtained for both normal and CML cells (figures 2(e) and (f)) (see appendix C). To this end, the p.d.f.s corresponding to ductile and brittle rupture events were disentangled with a two component Gaussian mixture model (see appendix B). When fixing ΔE₀ = 12 k_BT and ΔE* = 4 k_BT (characteristic unbinding energy of ABPs [46]), sets of parameters (a₀, $\hat{a}$, Δa) values (table 1) were found providing quite satisfactory fits of the experimental released energy distributions (figure 6). For ductile rupture events, the mean number of cascade steps is rather limited for normal ($\overline{N}=9$) as well as CML ($\overline{N}=8$) cells and correspond to low mean released energy values $\overline{E}$ and likely to a few tens of ABP unbindings (${\overline{N}}_{{\rm{unb}}}\sim \overline{E}/(4$ k_BT)) for normal ($\overline{E}=191$ k_BT, ${\overline{N}}_{{\rm{unb}}}\sim \,48$) and CML ($\overline{E}=198$ k_BT, ${\overline{N}}_{{\rm{unb}}}\sim \,50$) cells. For brittle failures, our model predicts more expanded cascades with significantly larger values of $\overline{N}$, $\overline{E}$ and ${\overline{N}}_{{\rm{unb}}}$ for both normal ($\overline{N}=21$, $\overline{E}=1104$ k_BT, ${\overline{N}}_{{\rm{unb}}}=276$) and CML ($\overline{N}=22$, $\overline{E}=1494$ k_BT, ${\overline{N}}_{{\rm{unb}}}\sim 373$) cells (table 1). Besides confirming that brittle rupture events involve more ABP unbindings than ductile ones, our rupture cascade model also predicts that the rate of energy released during the first cascade steps $\overline{{\rm{\Delta }}{E}_{n}}$ and the maximum reached before the exponential decrease to zero for large n are much higher for brittle than for ductile rupture events (figure 7(a)). Interestingly, when adimensionalizing ΔE_n by a mean released energy E/N per step and n by the total number of cascade steps (figure 7(b)), the numerical data obtained for ductile events in normal and CML cells superimpose and display a slow increase from ∼0.9 to ∼1.3 during the first part of the cascade before decaying toward the energy cut off. For brittle events, the numerical data for normal and CML cells again superimpose (figure 7(b)), but now exhibit an initial three-fold increase from ∼0.5 to ∼1.5. This initial acceleration of ABP unbindings from ∼2–3 to reach values as high as ∼15–20 unbindings per cascade step in brittle events confirms the existence of some correlation between the ABP unbindings of successive cascade steps that likely results in a collective local disorganization and possible disintegration of the CSK. In comparison, the smoother released energy cascade with only a few ABP unbindings (up to 5) at each cascade step strengthens our interpretation of ductile rupture events as more dynamical and reversible stress-induced cross-linker unbindings that would confer to the cell ductile plasticity to large deformations. We have simplified the discussion here by considering only one energy ≃4 k_BT for the actin cross-linker unbinding events, and grouping the brittle and ductile failures into two groups involving a different number of cross-linker unbinding events. Actually we could have also considered that the actin CSK network includes two populations of cross-linking proteins, namely weak cross-linking proteins that give a ductile failure and tight cross-linking proteins leading to brittle failures. The important physical quantity that is quantitatively predicted by the model is the total energy E released during a cascade of rupture events.

We refer the reader to the supplemental material (see footnote 11) where numerical simulations performed with different threshold ΔE₀ (12 and 16 k_BT) and arrest cut-off ΔE* (0 and 4 k_BT) (see tables S1–S3, and figures S9–S15 (see footnote 11)) confirm the robustness of our rupture cascade model predictions.

After informed consent in accordance with the Declaration of Helsinki and local ethics committee bylaws (from the Délégation à la recherche clinique des Hospices Civils de Lyon, Lyon, France), bone marrow samples were obtained from CP-CML patients at diagnosis and from healthy allogeneic donors. Mononuclear cells were separated using a Ficoll gradient (Bio-Whittaker) and were then subjected to CD34⁺ immunomagnetic separation (Stemcell Technologies). The purity of the CD34⁺ enriched fraction was checked by flow cytometry and was over 95% on average. Selected bulk CD34⁺ cells were seeded at 6 × 10⁵ cells ml⁻¹ and cultured in serum-free Iscove's Modified Dulbecco's Medium (Invitrogen) in the presence of 15% BSA, insulin and transferrin (Stemcell Technologies) supplemented with 10 ng ml⁻¹ interleukin-6 (IL-6), 50 ng ml⁻¹ stem cell factor, 10 ng ml⁻¹ IL-11 and 10 ng ml⁻¹ IL-3 (Peprotech).

To question the possible modifications of mechanical properties of immature hematopoietic cells upon oncogene expression, we have used the TF1 cell line transduced with the BCR-ABL oncogene as a unique model of CML cells that reproduces the early steps of stem cell transformation [40–44, 62–67]. In CML particularly, BCR-ABL is known for its ability to decrease progenitor cell adhesion to stroma [42, 68, 69]. By performing adhesion assays, we observed that BCR-ABL transduction has no effect on the percentage of adhesion of TF1 cells on fibronectin-coated surfaces (data not shown) [70]. This is according to the fact that TF1 cells are representative of very immature cells rather than progenitors of more differentiated cells. The immature CD34⁺ TF1 cell line (ATCC CRL-2003) was maintained at 1 × 10⁵ cells ml⁻¹ in RPMI-1640 medium, 10% fetal calf serum and granulocyte macrophage colony-stimulating factor (10 ng ml⁻¹) (Sandoz Pharmaceuticals). Engineered TF1-GFP and TF1-BCR-ABL-GFP cell lines were obtained by transduction with a murine stem cell virus-based retroviral vector encoding either the enhanced green fluorescent protein cDNA alone (EGFP) as a control or the BCR/ABL-cDNA upstream from an IRES-eGFP sequence [42, 64]. EGFP+ TF1 cells were sorted using a Becton Dickinson FACSAria.

Adhesion was performed as described in [68–70]. Fibronectin (Sigma) was coated for 2 h at 37 °C on 96-well Cellstar dishes (Greiner) using a 50 μg ml⁻¹ ligand solution in sodium bicarbonate buffer (0.1 mM). Plates were blocked with 0.3% BSA in PBS for 1 h at 37 °C. Cells (5.10⁴/well) were labeled for 20 min at 37 °C with 5 μM calcein-AM (Molecular Probes) in RMPI medium containing 0.1% BSA (without phenol red). Cells were then allowed to adhere to the coated plates in triplicate wells for 1 h at 37 °C in the adhesion buffer (PBS with 0.03% BSA), before removal of the non-adherent fraction. Finally, the glass coverslip was mounted on the AFM stage and the cells were kept in their culture medium at room temperature.

Cells were seeded at 5 × 10⁵ cells ml⁻¹ and incubated for 24 h before the experiment. On the one hand, non-adherent cells were centrifuged onto glass slides with a cytospin 4 (Thermo Scientific). On the other hand, non-adherent cells were allowed to adhere on glass cover slips coated with fibronectin (Sigma) in culture treated plates (BD Biosciences/Falcon) for 1 h at 37 °C, before removal of the non-adherent fraction. Cells were then fixed with 4% paraformaldehyde for 15 min, permeabilized with 0.1% Triton X-100 (Sigma-Aldrich) for 10 min, and blocked with 0.2% gelatin (Sigma-Aldrich) for 30 min at room temperature. The filamentous actin (F-actin) was labeled with 1000-fold diluted phalloidin-rhodamin (Sigma-Aldrich) for 30 min at room temperature. Cells were washed twice and were then incubated for 1 h at room temperature with a 400-fold diluted rabbit polyclonal anti-GFP antibody-Alexa Fluor 488 conjugate (Life Technologies). Finally, cells were washed twice and subsequently the nuclei were labeled with 2000-fold diluted DAPI (Sigma-Aldrich) in mounting medium (Sigma-Aldrich). Controls for nonspecificity and autofluorescence were performed by incubating cells without phallodin-rhodamin or by incubating TF1 cells with the rabbit polyclonal anti-GFP Alexa Fluor 488 conjugate. Fluorescence images were taken using a spectral confocal microscope TCS SP5 AOBS (Leica) (For more details see [45, 70]).

To investigate the structural transformation of the actin CSK consecutive to BCR-ABL transduction, we performed two different types of microscopy studies: (i) confocal microscopy [45, 70] on fixed TF1 and TF1-BCR-ABL cells where both F-actin and the nucleus were stained (figure 3), and (ii) quantitative phase microscopy [71–73] from which an important disorganization of the internal cell compartments was put into light by enhanced optical phase gradients. In BCR-ABL-transduced TF1 cells, juxtanuclear actin aggregates were found in almost 30% of the cells in addition to the cortical F-actin staining [45] (figure 3(b)). Very likely these structures were induced by the oncogene since they were rarely observed in the parental TF1 cell line (figure 3(a)). Consistently BCR-ABL-transduced cells showed fewer actin stress fibers as compared to TF1 cells in adhesion. The mechanical confinement of the TF1-BCR-ABL cells onto adherent fibronectin surfaces enlightens the interplay between the oncogene BCR-ABL and actin microfilaments, making their morphological changes quite visible by fluorescence microscopy (figure 3). These observations indicate that transducing TF1 cells with BCR-ABL likely disrupts their actin CSK cohesion and inhibit their crawling motility [45].

FICs were recorded from AFM cantilever deflection signals, when indenting the cell (approach curves) at fixed scan velocity V₀ = 1 μm s⁻¹ with a CellHesion system equipped with a 15−200 μm motorized stage, from JPK Instruments AG (www.jpk.com) mounted on an inverted microscope. We used pyramidal shape tip cantilevers (SNL-10, Bruker) with a nominal spring constant of 0.06 nN nm⁻¹. Prior to each experiment, the deflection sensitivity of the cantilever was estimated on fused silica and the cantilever spring constant k was calibrated by the thermal noise method in between 0.05 and 0.15 N m⁻¹ by directly recording their free fluctuations in buffer solution, computing their power spectrum distribution and fitting these curves with Lorentzian distributions [74]. Cells were prepared by letting suspended cells adhere on glass cover slips coated with fibronectin (Sigma) in culture treated plates (BD Biosciences/Falcon) for 1 h at 37 °C, and the non-adherent fraction was removed before mounting the coverslip on the AFM stage. During AFM recording, the cells were kept in their culture medium at room temperature (24 °C). To reduce variability, indentation was carried out at the center of each cell within 2 h after removing cells from the incubator, to probe the perinuclear CSK that is known to play a protective and mechanical confining role for the underlying nucleus and its multi-scale functions [75]. In all the experiments, cell sizes were evaluated from transmission microscopy images. To avoid these cells slipping away from the cantilever during the indentation, we immobilized them on a fibronectin-coated coverslip before AFM probing. Interaction of the integrins at the cell membrane with the fibronectin-coated surface not only confines cell movements on the glass surface but also modifies its CSK architecture, inducing a cascade of molecular events leading to cell spreading and higher adherence [76]. Actually, this fibronectin adhesion assay mimics to some extent the fact that, under physiological conditions, hematopoietic cells are not in suspension in their human bone marrow niche but actively adhere to the stroma and to proteins of the extracellular matrix [42, 68].

If we define Z as the distance of the cantilever tip to the sample surface: Z = Z_tip − Z_surf, the contact point Z_c corresponds to the tip position touching the sample surface without being deflected. Z_c is taken as origin of the FICs in figures 1(c) and (d). Overall, one approach/retract experiment lasts a few seconds, which is shorter than the characteristic remodeling time by active molecular motors. A key issue in the analysis of FICs is the determination of Z_c for very soft materials, like living cells. To master this practical signal-to-noise issue, we used a wavelet-based decomposition of the FICs and their derivatives [36, 39] that likely achieved some compromise between a too strong smoothing of the force derivative that would wipe out the non-contact to contact transition, and a too mild smoothing of the force that would suffer from a noise estimation of the contact points (see figures S2 and S3 (see footnote 11)). Once the cell is deformed by the cantilever, Z > Z_c, Z − Z_c reads as the sum of two terms, respectively the cell indentation h and the ratio of the force F and the cantilever spring constant k [77]:

$$\begin{eqnarray}&&Z-{Z}_{c}=h+F/k.\end{eqnarray} \tag{ A1 }$$

The nominal spring constant k of the cantilever (k = 0.06 nN nm⁻¹) was chosen large enough for the cantilever deflection to be small compared to the cell deformation. Nevertheless, we have substracted the correcting term from the FICs for energy integral computation. Note that this correcting term does not affect the computation of the second-order derivative of the force.

The wavelet transform is a mathematical time–frequency (time scale) decomposition of signals introduced in the early 1980s [78]. The wavelet transform has been applied to a great variety of situations in physics, physical chemistry, biology, signal and image processing, material engineering, mechanics, economics, epidemics... [75, 79–84]. Real experimental signals are very often nonstationary (they contain transient components), and when they are complex or singular, they may also involve a rather wide range of frequencies. It also happens that experimental signals display characteristic frequencies that drift in time. Standard Fourier analysis is therefore inadequate in these situations, since it provides only statistical information about the relative contributions of the frequencies involved in the analyzed signal. The possibility to perform simultaneously a temporal and frequency decomposition of a given signal was first proposed by Gabor for the theory of communication [85]. Later on, two distinct approaches (based on different wavelet transforms) were developed in parallel: (i) a continuous wavelet transform (CWT) [78, 82] and (ii) a discrete wavelet transform [80]. For singular (self-similar or multi-fractal) signals or images, the CWT transform rapidly became a predilection mathematical microscope to perform space (or time)-scale analysis and to characterize scale invariance properties. In particular it was used to generalize box-counting techniques [86] and to remedy the limitations of structure function methods [87] in elaborating a statistical physics formalism of multifractals [81, 84, 87–92].

During the past 30 years, the CWT was used for biological applications, on both 1D signals and 2D images [75, 83, 84, 93–107]. As far as 1D signals are concerned, the CWT was applied to AFM force curves collected from single living plant cells [39], living hematopoietic stem cells [36, 45] and to AFM fluctuation signals to characterize the passive microrheology of living myoblasts [108]. The 1D CWT was also generalized to 2D (and to 3D) CWT [79, 83, 84, 96] and it proved again its versatility and power for analyzing AFM topographic images of biosensors [109], fluorescence microscopy images of chromosome territories [98] and diffraction phase microscopy of living cells [71–73].

Viscoelastic theories developed during the second half of the twentieth century have led to general hereditary integral representation of stress-strain relationships for the indentation of linear viscoelastic materials by axisymmetric indenters [110]:

$$\begin{eqnarray}&&F(t)=\displaystyle \frac{4}{1-\nu }{C}_{n}{\int }_{0}^{t}G(t-\tau )\displaystyle \frac{{{\rm{d}}{h}}^{(n+1)/n}(\tau )}{{\rm{d}}\tau }{\rm{d}}\tau ,\end{eqnarray} \tag{ B1 }$$

where G(t) is the stress relaxation modulus, F is the loading force, h describes the displacement of the indenter and n is a positive integer which depends on the shape of the indenter. The stress relaxation modulus G(t) retains the strain-rate memory of the deformation. For a pyramidal indenter tip [38], we have n = 1 and ${C}_{1}=\tan \theta /\pi$, where θ is the nominal tip half-angle:

$$\begin{eqnarray}&&F(t)=\displaystyle \frac{4\tan \theta }{\pi (1-\nu )}{\int }_{0}^{t}G(t-\tau )\displaystyle \frac{{{\rm{d}}{h}}^{2}(\tau )}{{\rm{d}}\tau }{\rm{d}}\tau .\end{eqnarray} \tag{ B2 }$$

Since the cantilever is swept at constant velocity V₀, ${\rm{d}}{Z}={V}_{0}{\rm{d}}{t}$, the stress relaxation modulus G can be rewritten as:

$$\begin{eqnarray}&&G(Z)=\displaystyle \frac{\pi (1-\nu )}{8\tan \theta }\displaystyle \frac{{{\rm{d}}}^{2}F(Z)}{{{\rm{d}}{Z}}^{2}},\end{eqnarray} \tag{ B3 }$$

meaning that the change of G with Z is an important quantity that is stored in the whole deformation story. This relation establishes that the stress relaxation modulus can be obtained from the second derivative of the indentation force curve with respect to Z, without assuming a priori a particular viscoelastic or plastic cellular model.

Practically, we used a time–frequency adaptative wavelet-based method to compute G(Z) from FICs. Within the norm ${{ \mathcal L }}^{1}$, the one-dimensional WT of a signal F(Z) reads [78–82]:

$$\begin{eqnarray}&&{W}_{\psi }[F](b,s)=\displaystyle \frac{1}{s}{\int }_{-\infty }^{\infty }F(Z){\psi }^{* }\left(\displaystyle \frac{Z-b}{s}\right){\rm{d}}{Z},\end{eqnarray} \tag{ B4 }$$

where b is a spatial coordinate (homologous to Z) and s (> 0) a scale parameter. In the context of this study, we concentrated on the family of analyzing wavelets obtained from the successive derivatives of the Gaussian function [75, 81–84] (see figure S2(a) (see footnote 11)):

$$\begin{eqnarray}&&{g}^{(0)}(Z)={{\rm{e}}}^{-{Z}^{2}/2}.\end{eqnarray} \tag{ B5 }$$

Let us define the first derivative of the Gaussian function (see figure S2(b) (see footnote 11)):

$$\begin{eqnarray}&&{g}^{(1)}(Z)=-\displaystyle \frac{{\rm{d}}}{{\rm{d}}{Z}}{g}^{(0)}(Z)={{Z}{\rm{e}}}^{-{Z}^{2}/2},\end{eqnarray} \tag{ B6 }$$

and its second derivative, also called the Mexican hat wavelets (see figure S2(c) (see footnote 11)):

$$\begin{eqnarray}&&{g}^{(2)}(Z)=-\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{Z}}^{2}}{g}^{(0)}(Z)={{\rm{e}}}^{-{Z}^{2}/2}(1-{Z}^{2}).\end{eqnarray} \tag{ B7 }$$

Via one (resp. two) integration(s) by part, it is straightforward to demonstrate [75, 81–84, 90] that the WT of F with the first (resp. the second) derivative of a Gaussian wavelet at scale s, ${W}_{{g}^{(1)}}[F](b,s)$ (resp. ${W}_{{g}^{(2)}}[F](b,s)$) is precisely the first (resp. second) derivative of a smooth version ${W}_{{g}^{(0)}}[F](b,s)$ of F by a Gaussian function at the same scale s:

$$\begin{eqnarray}&&{W}_{{g}^{(1)}}[F](b,s)=s\displaystyle \frac{{\rm{d}}}{{\rm{d}}{b}}{W}_{{g}^{(0)}}[F](b,s),\end{eqnarray} \tag{ B8 }$$

and

$$\begin{eqnarray}&&{W}_{{g}^{(2)}}[F](b,s)={s}^{2}\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{b}}^{2}}{W}_{{g}^{(0)}}[F](b,s).\end{eqnarray} \tag{ B9 }$$

Let us point out that the validity of the WT definition (equation (B4)) was further extended for distributions including Dirac distributions [81, 82, 90, 111]. The interest of the WT method is two-fold. The first advantage is to use the same smoothing function to filter out the experimental background noise and to compute higher-order derivatives (for instance up to second-order in this study) at a well defined smoothing scale s_g. The second advantage relies on the power of the WT to detect local singularities (including rupture events in the FICs) and to quantify their force via the estimate of local ${\rm{H}}\ddot{o}$lder exponents from the behavior across scales of the WT modulus maxima [75, 81–84, 87–92]. In this study, we used modified versions [45] of the definition (equation (B4)) of the WT to get a direct measure of F in nN (${T}_{{g}^{(0)}}[F](b,s)$), dF/dZ in nN nm⁻¹ (${T}_{{g}^{(1)}}[F](b,s)$) and ${{\rm{d}}}^{2}F/{{\rm{d}}{Z}}^{2}$ in Pascal (${T}_{{g}^{(2)}}[F](b,s)$), once smoothed by a Gaussian window (${g}^{(0)}(Z)$) of width s:

$$\begin{eqnarray}&&{T}_{{g}^{(0)}}[F](b,s)={W}_{{g}^{(0)}}[F](b,s),\end{eqnarray} \tag{ B10 }$$

$$\begin{eqnarray}&&{T}_{{g}^{(1)}}[F](b,s)=\displaystyle \frac{1}{s}{W}_{{g}^{(1)}}[F](b,s),\end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray}&&{T}_{{g}^{(2)}}[F](b,s)=\displaystyle \frac{1}{{s}^{2}}{W}_{{g}^{(2)}}[F](b,s).\end{eqnarray} \tag{ B12 }$$

Figure S3 (see footnote 11) illustrates on a single FIC (approach and retract curves) the computation of the wavelet-based force derivatives. Let us point out that from equations (B1)–(B3), we get the following expression for the stress relaxation modulus [36, 39, 45]:

$$\begin{eqnarray}&&G(Z)=\displaystyle \frac{\pi (1-\nu )}{8\tan \theta }{T}_{{g}^{(2)}}[F](Z,s).\end{eqnarray} \tag{ B13 }$$

The protocol that we have elaborated to track singular events in FICs is shown in figures 1(e) and (f). The optimal wavelet scale to compute G(Z) via equation (B13) was fixed to $w=2\sqrt{2}s=42\,\mathrm{nm}$. Noticing that the FIC disruption events occurred in between two consecutive minima G_m and maxima G_M of ${{\rm{d}}}^{2}F(Z)/{{\rm{d}}{Z}}^{2}$ (figure 1(f)), we took the local minima G_m as searching criteria. We defined a threshold $| {G}_{m}|$ for the distribution of G_m values computed on a representative set of FICs, to discriminate the disruption events from the background noise. The prominence of these negative peaks was set to $| {G}_{m}| \geqslant 50$ kPa. In the right neighborhood of these peaks, we searched for a local maxima of ${{\rm{d}}}^{2}F(Z)/{{\rm{d}}{Z}}^{2}$ with a peak prominence $\geqslant 20$ kPa. G_m and G_M are marked by black symbols in figure 1(f). From the two positions of G_m and G_M we could then detect the beginning and the end of the disruption events represented by blue symbols in figure 1(e). The distance between these two positions was denoted as ΔZ. Finally, the force drop ΔF was corrected, taking into account the increasing trend of the FIC after the rupture event. A linear interpolation of the FIC in a small interval (∼30 nm) beyond the local minima of the FIC gave the best interpolation.

The normalized histograms obtained for ${\mathrm{log}}_{10}({\rm{\Delta }}Z)$, ${\mathrm{log}}_{10}({\rm{\Delta }}F)$ and ${\mathrm{log}}_{10}(E)$ for the rupture events detected in the two sets of FICs from healthy and CML cells were fitted (nonlinear least square fit method) by the sum of two Gaussian distributions:

$$\begin{eqnarray}&&P(Y)=\displaystyle \frac{(1\,-\,\alpha )\exp \,\left[\tfrac{-{(Y-{\mu }_{1})}^{2}}{2{\sigma }_{1}^{2}}\right]}{\sqrt{2\pi }{\sigma }_{1}}+\displaystyle \frac{\alpha \exp \,\left[\tfrac{-{(Y-{\mu }_{2})}^{2}}{2{\sigma }_{2}^{2}}\right]}{\sqrt{2\pi }{\sigma }_{2}},\end{eqnarray} \tag{ B14 }$$

where the indices 1 and 2 refer to ductile and brittle rupture events respectively. μ₁ (resp. μ₂) and σ₁ (resp. σ₂) are the mean and root-mean-square of the random variable $Y={\mathrm{log}}_{10}(X)$. The arithmetic and the geometry means of X are $\overline{X}=$ E $[X]={10}^{\mu +\mathrm{ln}(10){\sigma }^{2}/2}$ and $\tilde{X}=$ GM $[X]={10}^{\mu }$ respectively. The parameter α in equation (B14) quantifies the statistical contribution (%) of brittle events as compared to the contribution 1 − α (%) of ductile rupture events. In figure 6, a two component Gaussian mixture model was used to disentangle the p.d.fs of ${\mathrm{log}}_{10}(E)$ corresponding to brittle rupture events and ductile rupture events.