Computer simulations of the mechanical response of brushes on the surface of cancerous epithelial cells

We report a model for atomic force microscopy by means of computer simulations of molecular brushes on surfaces of biological interest such as normal and cancerous cervical epithelial cells. Our model predicts that the force needed to produce a given indentation on brushes that can move on the surface of the cell (called “liquid” brushes) is the same as that required for brushes whose ends are fixed on the cell’s surface (called “solid” brushes), as long as the tip of the microscope covers the entire area of the brush. Additionally, we find that cancerous cells are softer than normal ones, in agreement with various experiments. Moreover, soft brushes are found to display larger resistance to compression than stiff ones. This phenomenon is the consequence of the larger equilibrium length of the soft brushes and the cooperative association of solvent molecules trapped within the brushes, which leads to an increase in the osmotic pressure. Our results show that a careful characterization of the brushes on epithelial cells is indispensable when determining the mechanical response of cancerous cells.

The conservative force is given by a soft, linearly decaying repulsive function: where rij = ri − rj, rij = |rij|, ̂ = rij/rij, rij is the magnitude of the relative position between particles i and j, and is the intensity of the repulsion between a pair of particles. This interaction constant depends on the coarsegraining degree [S2], i. e., on the number of water molecules grouped into a DPD bead. In all our simulations we use a coarsegraining degree equal to three, which leads to = 78.0( ⁄ ) for all types of particles, except for the interaction between solvent particles and the monomers that make up the polymeric brushes; in such case, = 79.3( ⁄ ); is the Boltzmann's constant and T the absolute temperature. The dissipative and the random forces are, respectively: where is the noise amplitude and is the friction coefficient and they are related as follows: = 2 2 ⁄ ; vij = vi − vj is the relative velocity between the particles, and = is a random number uniformly distributed between 0 and 1 with Gaussian distribution and unit variance. The weight functions and depend on distance and vanish for r > , and are chosen for computational convenience to be (see [S1]): All forces between particles i and j vanish further than a finite cutoff radius , which represents the inherent length scale of the DPD model and it is regularly chosen as the reduced unit of length, = 1. The constants in equations (S3) and (S4) are chosen as  = 3 and  = 4.5, so that kBT = 1 [S3], for all the cases reported in this work because those constants fix the thermostat, which should be the same for all types of brushes modeled here.
The DPD particles that make up the brushes and the buffer in this work are confined by a simulation box which is flanked by effective square surfaces placed at the ends of the box, perpendicularly to the zaxis, i.e. at z = 0 and z = Lz. The force model for those surfaces is a simple, linearly decaying force law in the zdirection. It is given by: There are more sophisticated DPD wall models [S4], but for the purposes of the present work the model in equation (S6) is sufficient because it represents a soft wall, as a cell's surface is. It is important to emphasize that this model does allow for the deformation of cells' surfaces, as well as the brushes attached to them. The wall force in equation (S6)  also. As for the interaction of solvent and brush particles with the opposite surface, which is the tip of the AFM, they were all chosen as aw = 140.0 (kBT/ ). To model "liquid" brushes [S5] we allowed the attached ends of the polymeric brushes to move freely on the surface of the cell (the xy -plane), while for the modeling of "solid" brushes, those ends were fixed at random positions on the xyplane. The brush chains were modeled as linear chains made up of identical beads which are joined by freely rotating harmonic springs [S6]: where the spring constant, 0, was set at 0 = 100.0 (kBT/rc 2 ) for the so -called "soft" brushes, including those on normal cells; for the "stiff" brushes we chose 0 = 2000.0 (kBT/rc 2 ). The springs' position of equilibrium, r0, was chosen in all cases as r0 = 0.7rc [S7]. All beads that make up the polymer chains that model the brushes on normal and cancerous cells have the same size (rc), which is larger than the persistence length of the brush "molecules". Therefore torsional angles and vibrations of angles between bonds are important only within beads, and the relatively long chains that represent the microvilli become fully flexible and can be modeled accurately with equation (S7) [S8-S10].
The central purpose of this work is the calculation of the force that the tip of the AFM exerts on the surface brushes of normal and cancerous cervical cells, which is obtained through the disjoining pressure,  [S11]. It is defined as where (ℎ) is the component of the pressure tensor that is perpendicular to the xyplane, i. e., is the pressure exerted on the surfaces that confine the complex DPD fluid made up of polymer chains and solvent monomers. The bulk pressure, , is obtained from the average of the diagonal of the pressure tensor of the isotropic system, namely the one without confining surfaces. The components of the pressure tensor are calculated using the virial theorem route [S12], which provides kinetic and "potential" contributions to the pressure tensor. These components of the pressure tensor were calculated following the model of Irving and Kirkwood [S13], given as follows: where the first sum is the kinetic contribution. The second term is the product of the xcomponent of the conservative DPD force acting between particles i and j, and the xcomponent of the rij vector (see equation (S2)). The pressure tensor components Pyy and Pzz are obtained by replacing y and z by x in equation (S9), respectively. The normal pressure component in equation (S8)  as in some experiments [S14]. Because of the difference in scales of the brushes' thickness and the AFM probe one can use the socalled Derjaguin approximation [S15], where the force (F) per curvature radius (R) of the AFM tip can be obtained from the pressure in equation (S8) between flat surfaces as follows: Equation (S10) is valid when the distance between the two flat surfaces (h) is smaller than the curvature radius R, which is a condition satisfied by our simulations. The compression of the brushes on normal and cancerous cells was modeled in the range 2.8 nm ≤ h ≤ 17.5 nm and we used the value R = 2.5 m [S14], therefore the use of the Derjaguin approximation is fully adequate in this work. This is relevant for the scale of the experiments because these simulations focus on the scale where the "grafting density" of the liquid and solid brushes is equal, given that the tip of the AFM probes an area where there are brushes of different length and the force collected by the AFM is necessarily an average over such area.
We use a hybrid DPD-Monte Carlo (MC) algorithm in the grand canonical ensemble (fixed chemical potential, volume and temperature) as is required to reach the equilibrium state of all systems modeled under confinement [S11]. The dynamics part of the algorithm proceeds as follows: the motion of the DPD particles under the influence of forces given by equations (S1) - (S7) is solved through the integration of Newton's second law of motion [S12], using the velocity Verlet algorithm, adapted to the DPD interaction model [S16]. The dynamics is solved for ten time steps to generate a configuration, which is then compared with the initial configuration, after that the Metropolis algorithm is applied to choose between one of those two configurations [S17]. Once a configuration is chosen, the algorithm proceeds to try inserting and deleting aqueous solvent monomers a number of times, keeping the chemical potential of the solvent fixed. Full details of the DPD-MC method have been reported elsewhere [S11].

SYSTEMS STUDIED AND SIMULATION DETAILS
Dimensionless units are used throughout this work. Length is reduced with the cutoff radius, rc which, for a coarsegraining degree equal to three water molecules per DPD particle is equal to rc = 6.46 Å [S2]. The time step t is reduced with = ( 2 ⁄ ) 1 2 ⁄ , where m is the mass of a DPD particle; using the mass of three water molecules in a DPD bead, t = 6.3 ps at T=300 K. The reduced time step used in the dynamics part of the algorithm was in all cases t/(6.3 ps) = 0.02; the energy was reduced with kBT, with T = 300 K. The area of the simulation box was fixed in all cases with Lx/rc= 7, Ly/rc = 7, while Lz/rc was increased from 4.3 up to 27; periodic boundary conditions were used, except in the z-direction since this is the direction of the confinement. The chemical potential was fixed at µ=37.7kBT, which gives a total average numerical density of <>~3/rc 3 . The simulation results were obtained from averages after at least 600 blocks of 10 4 MC configurations each were carried out, with the first 100 blocks used to equilibrate the system and the rest were used for the production phase.
In each given block made up of 10 4 MC configurations, the percentage of successful MC moves was around (35 ± 2) %.
The brushes on normal cells were modeled as linear chains made up of the same number of beads, N = 27, with one of their ends "grafted" to the cell surface but free to move on it (the so -called "liquid" type of brush). The spring constant joining those beads was chosen for all these chains as 0 = 100.0 (kBT/rc 2 ), which is called a "soft" brush. There were 16 chains attached to the surface of the normal cell, which yields a grafting density of = 0.78 nm -2 .
Brushes on cancerous cells were designed as nonuniform, following the experimental clues [S14]; they were made up of three different types of brushes. The shortest one was made of 36 chains with N1 = 5 beads and grafting density 1 = 1.76 nm -2 ; the mediumsized brush consisted of 10 chains with N2 = 30 beads and 2 = 0.49 nm -2 , and the largest brush was made of 4 chains with N3 = 42 beads with 3 = 0.20 nm -2 . In contrast with the brushes on normal cells, the brushes on cancerous cells were allowed to have either fixed or moving ends on the surface (to model "solid" and "liquid" brushes, respectively), and the spring constant joining the beads along the chains was allowed to have two values: 0 = 100.0 (kBT/rc 2 ) for soft brushes, and 0 = 2000.0 (kBT/rc 2 ) for stiff ones.

SUPPLEMENTARY RESULTS AND DISCUSSION
In Fig. S1 we show the force that the AFM tip exerts on the brush that covers normal cells predicted by our simulations (circles in Fig. S1(a)) and those measured by Iyer and collaborators [S14] (circles in Fig. S1(b)). The trend predicted by our calculations is in very good agreement with that found in the experiments, although the scale of the distance between the AFM tip and the brush is orders of magnitude smaller. This is simply due to the  Force measured on normal human epithelial cervical cells (blue empty circles) with an AFM [S14]. The solid black line in both curves represents a best fit to the Alexanderde Gennes scaling equation for the force between polymer brushes [S15], see text for details of the fit.
Following Iyer and coworkers [S14] we have included in Fig. S1 the best fit to the Alexanderde Gennes scaling law for the osmotic pressure between a brush of uniform length and a solid surface [S18], which is given by: where R is the radius of curvature of the AFM tip,  is the polymers' grafting density and L is the equilibrium length of the brush, namely its length when it is not compressed by the AFM tip. Equation (S11) is expected to be valid only when 0.1 < ℎ < 0.9 ⁄ . We have used equation (S11) to fit our data and the result is the black line in Fig. S1(a), which is given by the equation (ℎ) = (540 nN) −2 ℎ⁄ with an effective equilibrium length L = 8 nm; in reduced DPD units this length is L/rc= 80/6.46 ≈ 12.4, which is in good agreement with the equilibrium length obtained for the brushes in normal cells from simulations, as is indicated in Fig. S2. Our simulations have the advantage that the chains that make up the brushes do interact with one another, and do not assume that the brush density profile (Fig. S2) is a step function [S18]; none of these features are incorporated into the approximate scaling function in equation (S11) [S18]. Regarding our best fit of the data of Iyer and collaborators [S14] using equation (S11) (black line in Fig. S1(b)), we obtain (ℎ) = (24.4 nN) −2 ℎ⁄ , with equilibrium length L = 2175 nm, which is to be compared with the average experimental value, see the Supplementary Information in reference [S14]: 2360 ± 970 nm. Therefore, DPD simulations are a powerful tool to study the mechanical response of brushes on cells. Figure S2. Density profile of a brush on a normal epithelial cell predicted by simulation. In this figure we see the number density of the beads that make up the chains of uniform length (solid squares), which is our model for the brush on normal cells, when the tip of the AFM is far from the brush, so that it can be considered to be in equilibrium. The equilibrium length of the brush (L) is indicated by the red line. Both axes are shown in reduced DPD units.
In Fig. S3 we show the density profiles of the soft and stiff liquid brushes in equilibrium, i.
e., when the tip of the AFM is far from the brushes so that their length is not perturbed by the confining force of the tip. As we have argued in the main manuscript, in the stiff brush the beads that form the chains are closer to each other than in the soft brush, because the spring is "harder", hence the equilibrium length of the stiff brush (Lstiff) is shorter than that of the soft brush (Lsoft). The soft brush then "swells" more than the stiff brush when molecules of the aqueous solvent penetrate it, leading to a larger force required to deform the soft brush 11 with the AFM than that required to deform the stiff brush. This argument follows also from equation (S11).

Figure S3. Comparison of brush density profiles on soft and stiff cancerous cells.
Number density profiles of the beads that make up the brushes of nonuniform density, joined by stiff and soft harmonic springs, when the brushes are in equilibrium, that is, when the tip of the AFM is not in contact with any chain in the brushes. The ends of the chains of the brushes that are attached to the cell surface are allowed to move on the plane of the surface, i. e., the brushes are liquid. Lsoft and Lstiff represent the equilibrium length of the soft and stiff liquid brushes, respectively.
In the experiments performed by Iyer and coworkers [S14] it was observed that cancerous cervical epithelial cells appeared to be covered by two rather dense brushes of different length, and a third brush of larger length than the other two, but of smaller "grafting density".
Their force profiles also seemed to indicate the presence of such third brush, although the authors fit their data using only a twobrush Alexanderde Gennes type of model, simply adding two terms of the type shown in equation (S11) with different parameters, see [S14].
In Fig. S4 we have followed the same approach to interpret the force profiles obtained from our simulations which, as we shall discuss, leads to important physical insights. h/r c Figure S4. Force profiles of brushes on soft and stiff solid brushes on model cancerous cells. In this figure we show only the force per curvature radius of the AFM tip obtained for the socalled solid brushes on cancerous cells, i. e., those brush "molecules" that have fixed attaching sites on the surface of the cell. The case corresponding to the liquid brushes was discussed in the previous figure and in the main manuscript. Filled, green triangles are the data obtained for soft brushes, while filled (red) circles represent the force per radius on stiff brushes. The dashed lines are the best fits to a three brush, Alexanderde Gennes type of model (see equation (S11)). See text for details.
In Fig. S4 we show the force profiles corresponding to the socalled solid type of brush only, since the characteristics of the liquid type have been discussed in detail in the previous paragraphs. The green triangles in Fig. S4 correspond to the force per curvature radius of the AFM tip when it compresses a soft, solid brush, while the red circles represent the equivalent profile, corresponding to the stiff solid brush. The dashed lines in Fig. S4 represent the best fits to a threebrush, Alexanderde Gennes type of scaling for the force per curvature radius, namely: In equation (S12), Ai and Li (with i = 1, 2, 3), x0, and B are adjustable parameters. After performing regression analysis, we obtained the following values of these parameters for the = 0.287 ± 0.003; A3 = 3.91, L3 = 7.0 ± 0.5, with x0 = 4.23 and B = -1.2 ± 0.1. As can be observed in Fig. S4, the fit of equation (S12) to our predictions is excellent. Notice that the values obtained for the maximum compression of the brushes from this fit (x0 = 4.28, and x0 = 4.23 for soft and stiff brushes, respectively) are in very good agreement with the actual value from our simulations (x0 = 4.3). More importantly, the fit yields an equilibrium length for the soft brush (Lsoft = L3 = 9.8 ± 0.4) that is larger than the equilibrium length of the stiff brush (Lstiff = L3 = 7.0 ± 0.5), which confirms our conclusion that the force required to compress a soft brush is larger than that required to compress the same distance a stiff brush.
We have shown here that this conclusion does not depend on the liquid or solid nature of the brushes, but only on the equilibrium length of the largest brushes the cover cancerous epithelial cells. Also, this work represents a successful coarsegrained model of an AFM cantilever interacting with a complex brush, which can be particularly useful when interpreting experiments on matter at the nanoscale, a task that still remains most challenging [S19].