Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modeling biomechanical responses of cells to external forces Wu, Tenghu 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2015_september_wu_tenghu.pdf [ 5.16MB ]
JSON: 24-1.0166427.json
JSON-LD: 24-1.0166427-ld.json
RDF/XML (Pretty): 24-1.0166427-rdf.xml
RDF/JSON: 24-1.0166427-rdf.json
Turtle: 24-1.0166427-turtle.txt
N-Triples: 24-1.0166427-rdf-ntriples.txt
Original Record: 24-1.0166427-source.json
Full Text

Full Text

Modeling biomechanical responses ofcells to external forcesbyTenghu WuB.Sc., Zhejiang University, 2008M.Sc., Zhejiang University, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Chemical and Biological Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)July 2015c© Tenghu Wu 2015AbstractFor many cells, their biomechanical properties are important to theirbiofunctions. This thesis contains three computational studies of cellulardynamics under mechanical deformation.When infected by malaria, infected red blood cells (iRBCs) become lessdeformable and tend to block microcapillaries. Microfluidic channels havebeen used to investigate the deformability of iRBC at different infectionstages. In my first project, I applied a discrete iRBC model to simulatethe traverse of iRBCs through a microfluidic channel and investigated theprogressive loss of the cell deformability due to three factors: the membranestiffening, the cell surface-volume ratio reduction, and the parasite growinginside the cell. The results indicate that the growth of the parasite clustersplay the most significant role in causing the channel blockage.Recent experiments have investigated the response of neutrophils af-ter passing through microfluidic channels. The results indicate that neu-trophils may be activated by mechanical deformation. Mechanical deforma-tion causes disassembly of the cytoskeletal network of the neutrophils, whichresults in a sudden drop of the cell elastic modulus (termed fluidization).The fluidization is followed by either activation of the neutrophils with for-mation of pseudopods or uniform recovery of the cytoskeletal network with-out pseudopod formation. The former only occurs when the neutrophils’transit rate is slow. I proposed a chemo-mechanical model for the fluidiza-tion and activation processes, based on the polarization of the Rac proteinthrough a wave-pinning mechanism. The model captures the main featuresof the experimental observation.The third project investigates the response of smooth muscle cells totransient stretch-compress (SC) and compress-stretch (CS) maneuvers. Priorexperimental results indicate that the transient SC maneuver causes a sud-den fluidization of the cell while the CS maneuver does not. To under-stand this asymmetric behavior, I built a biomechanical model to probethe response of stress fibers to the two maneuvers. The model couples thecross-bridge cycle of myosin motors with a viscoelastic Kelvin-Voigt element.Simulation results point to the sensitivity of the myosin detachment rate toiiAbstracttension as the cause for the asymmetric response of the stress fiber to theCS and SC maneuvers.iiiPrefaceThis thesis entitled “Modeling biomechanical responses of cells to ex-ternal forces” presents the research that I performed during my PhD study.The research conducted in this thesis was identified, initiated and supervisedby Professor James J. Feng. In this preface, the contributions and collab-orations to the papers published or submitted for publication from currentthesis are briefly explained.• A version of chapter 3 has been published: T. Wu and J. J. Feng(2013), Simulation of malaria-infected red blood cells in microfluidicchannels: Passage and blockage. Biomicrofluidics 7, 044115. In thisproject, I developed the numerical model, performed simulations andcollected data. The analyses was done by me and Professor James J.Feng together.• A version of chapter 4 has been submitted for publication: T. Wu andJ. J Feng (2015), Modeling the mechanosensitivity of neutrophils pass-ing through microfluidic channels. I developed the neutrophil model,performed all simulations and collected data. I and Professor JamesJ. Feng drafted the manuscript for this project.• A version of chapter 5 has been published: T. Wu and J. J. Feng(2015), A biomechanical model for fluidization of cells under dynamicstrain. Biophysical Journal 108, 43–52. I developed the mathematicalmodel for the stress fibers, performed all simulations and collecteddata. I and Professor James J. Feng analyzed the results and wrotethe paper together.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and objectives of this thesis . . . . . . . . . . . . 11.2 Mechanical properties of cells . . . . . . . . . . . . . . . . . . 21.3 Cellular mechanosensing . . . . . . . . . . . . . . . . . . . . 31.4 Experimental techniques . . . . . . . . . . . . . . . . . . . . 41.4.1 Micropipette aspiration . . . . . . . . . . . . . . . . . 41.4.2 Microfluidic channels . . . . . . . . . . . . . . . . . . 41.4.3 Substrate stretcher . . . . . . . . . . . . . . . . . . . 51.5 Cell models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.5.1 Mechanical models . . . . . . . . . . . . . . . . . . . 51.5.2 Chemo-mechanical models . . . . . . . . . . . . . . . 82 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Smoothed Particle Hydrodynamics (SPH) . . . . . . . . . . . 102.1.1 Kernal integral representation . . . . . . . . . . . . . 102.1.2 Particle approximation . . . . . . . . . . . . . . . . . 122.1.3 Kernel smoothing function . . . . . . . . . . . . . . . 132.1.4 Weakly compressible SPH algorithm . . . . . . . . . . 142.1.5 Boundary treatment . . . . . . . . . . . . . . . . . . . 16vTable of Contents2.2 Immersed Boundary Method . . . . . . . . . . . . . . . . . . 162.2.1 Governing equations . . . . . . . . . . . . . . . . . . . 162.2.2 Imposition of immersed boundaries . . . . . . . . . . 172.2.3 The smoothed δh function . . . . . . . . . . . . . . . 182.2.4 Computational algorithm . . . . . . . . . . . . . . . . 192.3 Cell models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3.1 Red blood cell (RBC) model . . . . . . . . . . . . . . 202.3.2 Neutrophil cell model . . . . . . . . . . . . . . . . . . 202.3.3 Stress fiber model for adherent cells . . . . . . . . . . 203 Simulation of malaria-infected RBC in microfluidic chan-nels: Passage and blockage . . . . . . . . . . . . . . . . . . . . 213.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Physical model and numerical scheme . . . . . . . . . . . . . 233.2.1 Cell membrane elasticity . . . . . . . . . . . . . . . . 243.2.2 Parameters for the model . . . . . . . . . . . . . . . . 253.2.3 Cell geometries . . . . . . . . . . . . . . . . . . . . . 263.2.4 Repulsive force between membrane and wall . . . . . 293.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.3.1 Effect of the membrane stiffness . . . . . . . . . . . . 323.3.2 Effect of the cell geometry . . . . . . . . . . . . . . . 353.3.3 Effect of the parasites cluster . . . . . . . . . . . . . . 383.3.4 Combinations . . . . . . . . . . . . . . . . . . . . . . 393.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414 Modelling the mechanosensitivity of neutrophils passing throughmicrofluidic channels . . . . . . . . . . . . . . . . . . . . . . . . 434.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.2 Neutrophils model . . . . . . . . . . . . . . . . . . . . . . . . 464.2.1 Mechanics of membrane deformation . . . . . . . . . 464.2.2 Kinetics of chemical signaling . . . . . . . . . . . . . 494.2.3 Remodeling of the cytoskeleton . . . . . . . . . . . . 514.3 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . 524.3.1 Fluid dynamics . . . . . . . . . . . . . . . . . . . . . 524.3.2 Immersed boundaries . . . . . . . . . . . . . . . . . . 524.4 Simulation setup and model parameters . . . . . . . . . . . . 534.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 554.5.1 Fluidization of the neutrophils . . . . . . . . . . . . . 564.5.2 Polarization of the neutrophils . . . . . . . . . . . . . 574.5.3 Activation of the neutrophils . . . . . . . . . . . . . . 66viTable of Contents4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 715 Fluidization of cells under dynamic strain . . . . . . . . . . 735.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.2 Model for stress fibers . . . . . . . . . . . . . . . . . . . . . . 755.3 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . 795.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 825.4.1 Responses to the SC and CS maneuvers . . . . . . . . 825.4.2 Stretch-hold (SH) maneuver . . . . . . . . . . . . . . 865.4.3 Critical conditions . . . . . . . . . . . . . . . . . . . . 895.4.4 Comparison with experiments . . . . . . . . . . . . . 925.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 936 Conclusions and recommendations . . . . . . . . . . . . . . . 966.1 Malaria-infected red blood cells . . . . . . . . . . . . . . . . 966.2 Activation of neutrophils . . . . . . . . . . . . . . . . . . . . 976.3 Cell fluidization and reinforcement . . . . . . . . . . . . . . . 97Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99viiList of Tables3.1 Dimensions of the healthy and infected red cell in differentstages. se = S/(4pia2) denotes the excess surface area ratioof the cell, where a is the effective cell radius defined as theradius of a sphere having the volume V of the cell. The lastcolumn indicates the dimensions of the parasite. . . . . . . . 263.2 The passage and blockage of healthy and infected red bloodcells through microfluidic channels of different size. . . . . . . 404.1 Parameters for the model. The value for the kpoly is calibratedthrough comparing the simulation with the experiment (seediscussion in Subsection 4.5.3). . . . . . . . . . . . . . . . . . 545.1 Parameters for the KVM stress-fiber model. . . . . . . . . . . 81viiiList of Figures2.1 The support domain of a particle at r. . . . . . . . . . . . . . 112.2 (a) Diagram for a solid body immersed in a flow. The bodytakes the volume Ωb with boundary Γb. (b) Schematic of bodyimmersed in a Cartesian grid. From Mittal and Laccarino[138] with permission, c©2005, Annual Reviews. . . . . . . . . 172.3 (a) Distribution of forcing F k from Lagrangian boundarypoint (Xk) to surrounding fluid nodes. Shaded region de-notes the force distribution region. (b) Distribution functionsemployed in various studies. From Mittal and Laccarino [138]with permission, c©2005, Annual Reviews. . . . . . . . . . . . 193.1 Passage and blockage of microfluidic channels of different sizesby iRBC in the ring, trophozoite and schizont stages of infec-tion. The flow goes from the right to the left. From Shelbyet al. [195] with permission, c©2003 The National Academyof Sciences. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 The healthy RBC and the ring-stage iRBC have the samebiconcave shape described by Eq. (3.6), with the smallest andlargest thickness being T1 = 0.81 µm and T2 = 2.4 µm. Thespatial coordinates are marked in microns. . . . . . . . . . . . 273.3 The cell geometries at (a) the early trophozoite stage, with(D1, D2, T1, T2) = (7.77, 7.76, 2.7, 3.0) µm, and (b) thelate trophozoite stage, with (D1, D2, T1, T2) = (7.80, 7.72,2.8, 4.0) µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4 Close up of Hemidactylus sp. . . . . . . . . . . . . . . . . . . 293.5 Schematic of the computational domain. The overall dimen-sions of the channel are L = 24 µm, W = 12 µm andH = 4.8 µm, with three segments of equal length l = 8 µm.The narrow section in the middle has a width w = 3.2, 4 and4.8 µm for the narrow, medium and wide channels. . . . . . . 30ixList of Figures3.6 The trajectories of RBCs through the narrow channel (w =3.2 µm) for three values of the shear modulus Gs. The insetsshow top-view snapshots of the cell during the transit withcolor contours of the local stretching γ of Eq. (3.8). Thered arrow pointed insets are for Gs = 15 µN/m at t = 0.27and t = 1.9, while the blue arrow pointed inset is for Gs =25 µN/m at t = 2.1. . . . . . . . . . . . . . . . . . . . . . . . 333.7 The transit time tt, made dimensionless by T , as a functionof the membrane shear modulus Gs for the three channels. . . 343.8 Steady-state shape of the iRBC in the early trophozoite, latetrophozoite and schizont stages after it blocks the narrowchannel (w = 3.2 µm). The length of the tongue lt, definedas the distance between the cell front and the channel entry,and scaled by the channel length L, is plotted as a function ofthe excess surface area ratio se. Color contours of the surfacestretching γ are also shown in the insets. . . . . . . . . . . . . 363.9 The transit time tt through the medium and wide channels asfunctions of the cell excess surface area ratio se. The transittime is scaled by the T . . . . . . . . . . . . . . . . . . . . . . 373.10 Trajectories of the schizont-stage iRBC in the wide channel(w = 4.8 µm) with a parasite of different sizes. R = 0 indi-cates a baseline case without parasite. . . . . . . . . . . . . . 394.1 Microfluidic channel. Cited from Yap and Kamm[238]. . . . . 454.2 (a) Schematic of the particle-spring model for the cell mem-brane. The portion inside the red box is shown in the magni-fied view on the right. (b) Xi denotes the location of the ithparticle, li is the length of the ith edge, and Θi is the anglebetween the ith and i+1th edges. At the start of simulations,all edges have the same length l0, and all neighbouring edgeshave the same angle Θ0 = 2pi/N . . . . . . . . . . . . . . . . . 474.3 Geometry for the channel with a constricting section. . . . . . 534.4 Passage of the neutrophil through the channel for γ = 4 s−1at (a) t = 1.54 s, (b) t = 2.30 s, (c) t = 3.07 s, (d) t = 3.84s, (e) t = 4.61 s and (f) t = 5.38 s. The colors represent thecortex modulus. . . . . . . . . . . . . . . . . . . . . . . . . . . 574.5 Passage of the neutrophil through the channel for γ = 4 s−1at (a) t = 1.54 s, (b) t = 2.30 s, (c) t = 3.07 s, (d) t = 3.84s, (e) t = 4.61 s and (f) t = 33.28 s. The colors indicate thelevel of active Rac a on the membrane. . . . . . . . . . . . . . 59xList of Figures4.6 Evolution of the total activated signal for γ = 2 s−1 (solidline) and γ = 4 s−1 (dash line). . . . . . . . . . . . . . . . . 604.7 Evolution of the signal at the cell’s front tip (af , solid line)and back tip (ab, dashed line) for (a) γ = 2 s−1 and (b) γ = 4s−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.8 The null-cline for the simplified reaction-diffusion equation(Eq. 4.19). Trajectory (i) corresponds to the polarizing case(γ = 4 s−1) while trajectory (ii) to the non-polarizing case(γ = 2 s−1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.9 Evolution of the signal content at the cell front (solid line)and back (dash line) for γ = 3 s−1. . . . . . . . . . . . . . . . 644.10 Phase diagram for three regimes of polarization, in the pa-rameter space of (a) ∆P versus te, and (b) ∆P versus γte. . 654.11 Distribution of F-actin inside activated neutrophils after pass-ing through filters. The cells are stained with tetramethyl-rhodamine isothiocyanate-phalloidin, and exhibit highly lo-calized concentrations of F-actin at sites of pseudopod protru-sion. (a) shows a neutrophil filtered through 3-µm pores un-der a constant pressure drop, while (b) show a neutrophil fil-tered through 3-µm pores under a constant flow rate. Adaptedfrom Yap and Kamm [237] with permission; c©the AmericanPhysiological Society. . . . . . . . . . . . . . . . . . . . . . . . 674.12 Evolution of the cell deformation index (DI) for the slow ki-netic rate γ = 2 s−1 and the fast rate γ = 4 s−1. The snap-shots (a)-(d) show the cell shapes at four interesting times. . 684.13 Recovery of the cortex elasticity for (a) the inactive case atslow kinetic rate γ = 2 s−1, and (b) the activated case atfaster kinetic rate γ = 4 s−1. In each case, we plot the tem-poral evolution of the cortical modulus at the front and reartips of the cell and the average over the entire membrane. Allthe moduli are normalized by the resting state modulus E0c . . 694.14 Predicted geometry and polarity of the activated neutrophilfor γ = 4 s−1 at t = 33.3 s. (a) Distribution of the corticalmodulus on the membrane. (b) Distributions of the Rac con-centration a and the protrusion force Fpro indicated by thearrows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70xiList of Figures5.1 Our model for the stress fiber has a Kelvin-Voigt elementconnected in series to a myosin apparatus. The two ends ofthe stress fiber are attached to the substrate through focaladhesions, which are assumed permanent in our model. . . . . 765.2 (a) The swinging crossbridge model adapted from Spudlich[201]. (b) A simplified three-stage model represents the in-teraction between the myosin and actin filaments, with m0,m1 and m2 denoting myosin in the detached, attached pre-stroke and attached post-stroke states. These states corre-spond roughly to those marked in (a). . . . . . . . . . . . . 775.3 Tension in the stress fiber under the SC (a) and CS (b) ma-neuvers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 825.4 Evolution of myosin populations under the SC maneuver. (a)m1 and m2; (b) the total attached myosin ma = m1 +m2. . . 845.5 Evolution of myosin populations under the CS maneuver: (a)m1 and m2; (b) the total attached myosin ma = m1 +m2. . . 855.6 Evolution of the tension σ during an SH maneuver modeledafter the experiment of Trepat et al. [213]: a 10% strain isapplied over T/2 = 2 s and then maintained indefinitely. Thecharacteristic time T = 4 s. ∆σ is the excess tension aftercessation of the stretch, and td is the time for ∆σ to relax by50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.7 Comparison of the excess tension, upon cessation of stretch,between the model prediction and experimental data. . . . . 885.8 Critical conditions for SF disassembly depicted on the (ε0, ε˙0)plane. For the experimental data, open symbols correspondto fluidization, filled ones non-fluidization, and half-filled onesthe critical condition. (a) The SC maneuver; (b) The CS andCH maneuvers. Model prediction at a lower Vm = 0.025 µm/sis also shown for comparison. . . . . . . . . . . . . . . . . . . 91xiiAcknowledgementsI would like to appreciate my supervisor, Professor James J. Feng. Thisdissertation could not be finished without his support and supervision. Ilearnt a lot from him during my PhD study. He taught me how to becritical on my research, and encouraged me when I got disappointed frommy work. I also would like to thank Professor Leah Edelstein-Keshet forgiving me many valuable advices during my PhD studies.xiiiDedicationThis thesis is dedicated to my parents. Thanks for their selfless love andsupport. I also want to express my sincere gratitude to my lovely girlfriend,Huimin Yun.xivChapter 1Introduction1.1 Motivation and objectives of this thesisIn my PhD thesis, I will investigate the biomechanical response of bio-logical cells to mechanical forces. It consists of three research projects, eachof which is inspired by recent experimental studies.• Shelby et al. [195] constructed a microfluidic channel to test the de-formability of malaria-infected red blood cells (iRBCs) at differentinfection stages. Their results imply a progressive stiffening of theiRBCs during the infection, as a result of modifications of the cell sub-structures. Such a device is promising for detecting malaria-infectedcells. The objective of this project is to build a mechanical model tosimulate the iRBC traversing a microfluidic channel, and to see howstructural modifications due to malaria lead to the rigidification of theiRBCs and blockage of the microfluidic channel.• Yap and Kamm [238] investigated neutrophil activation due to me-chanical deformations by pushing the cell through a microfluidic chan-nel. They found that once the neutrophil enters into the channel itselastic modulus initially drops (called fluidization) due to disassemblyof the cytoskeletal network. This is followed by formation of pseu-dopods due to polymerization of the actin filaments (called activa-tion). The fluidization seems to be due to mechancial perturbations.However, the activation is regulated through biochemical signals. Yapand Kamm [237] further found out that the activation also dependson the rate at which the cell enters the channel. The objective of thisproject is to build a chemomechanical model to simulate this inter-esting process, and to understand how such mechanical perturbationslead to fluidization and subsequently activation.• Trepat et al. [213] and Krishnan et al. [109] investigated how smoothmuscle cells respond to a transient stretch-compress deformation. Inthe experiment, the cells were cultured on an elastic substrate and the11.2. Mechanical properties of cellssubstrate was stretched by 10% strain and then released. The wholemaneuver lasts 4 s. They found that the maneuver causes a suddendrop of the elastic modulus of the cells (also called fluidization) dueto disassembly of the cytoskeleton. The fluidization happens at thebegin of the compressing phase. After the fluidization, the elasticmodulus of the cell gradually recovers to the resting level. Moreover,Chen et al. [18] found that a compress-stretch maneuver with 10%strain magnitude did not cause cell fluidization. The objective of thisproject is to build a biomechanical model to simulate this process, andto give a biophysical explanation for the distinct responses of the cellto the stretch-compress and compress-stretch maneuvers.1.2 Mechanical properties of cellsCells are the basic functional units of living organisms. Their mechanicalproperties are important to the biological functions of many cells [205, 206].Some diseases are induced by modifications of the cell’s mechanical proper-ties. Take malaria for example. Healthy red blood cells are extremely flexi-ble, which allows them to pass through micro-capillaries much smaller thantheir size. After infection by the malaria parasite Plasmodium falciparum,the cell gradually loses its deformability with the progress of infection due toparasitogenic proteins stiffening the cell membrane, growth of the parasitecluster, and reduction of the cell’s surface area to volume ratio. This causesblockage of capillary flows and the major symptoms of malaria.Generally speaking, the cell’s mechanical properties are determined by itscomponents, which include the plasma membrane, cytoskeleton, cell nucleusand the cytosol.• The plasma membrane is a lipid bilayer that separates the interior ofthe cell from its surrounding fluid. The plasma membrane is extremelyresistant to dilation deformation but does not resist shear deformation[141]. The plasma membrane also possesses bending force due to asym-metric lipid molecular density between the inner and outer lipid layers[72].• The cytoskeleton comprises actin filaments (F-actin) that are con-nected through actin-actin-binding proteins [232] and stressed by myosinmotors [175]. The cytoskeleton is usually attached to the plasma mem-brane through actin-membrane-binding proteins [232]. Under different21.3. Cellular mechanosensingexternal forcing conditions, the cytoskeleton may exhibit varied archi-tectures [57]. In turn, different cytoskeletal structures lead to differentmechanical properties of the cells. In suspended cells, the cytoskele-ton usually assumes the form of a 3D filamentous matrix, which iscommonly called the cortical network [148]. For adherent cells, thecytoskeleton usually takes on the form of contractile stress fibers (SFs)when the substrate is stiff [37]. The ends of the SFs connect to focaladhesions to sense external forces. In response to external stimuli, thecytoskeleton may undergo dramatic morphological changes that areregulated through the Rho GTPases [70].• The nucleus is the largest organelle in the cell. The average size of thenucleus for mammal cells is about 6 µm, which occupies about 10%of the total cell volume [2]. Typically, the nucleus appears in a dense,roughly spherical shape. The nucleus not only regulates the cell’sbio-functions, but also significantly influences the cell’s mechanicalproperties.• The cytosol is a viscous solution within the cell, mainly composed ofwater, salts, and proteins. Many other organelles are suspended in thecytosol, of course. But for our purpose we absorb their mechanicaleffects into the effective viscosity of the cytosol, and will not accountfor them explicitly.1.3 Cellular mechanosensingRecently, more and more evidence has revealed that the cell can respondto external mechanical forces with modification of its cytoskeletal structure,which in turn changes cell mechanical properties [109, 213, 238]. This phe-nomenon is termed cellular mechanosensing [122].Generally speaking, the mechanosensing process can be divided intothree steps: mechanical forces acting on the cell, mechanosensor proteinstransfer mechanical signals into biochemical signals, and biochemical sig-nals regulate the cytoskeletal dynamics and the cell responses. In vivo, cellslive in extremely complex physical and mechanical environments. Whenexternal forces initially act on the cell plasma membrane, they are trans-mitted to the cytoskeleton through the actin-membrane-binding proteins.The cytoskeleton is the major structure that balances the external forces.With the mechanical forces/strains acting on the cytoskeleton, mechanosen-sor proteins transform the mechanical signals into biochemical signals. The31.4. Experimental techniquesmechanosensor proteins typically consist of actin-membrane-binding [82] andactin-actin-binding [122] proteins. On the one hand, these proteins bindactin filaments into the cytoskeletal network and connects the cytoskele-ton to the plasma membrane. On the other hand, they can act as sensorsfor the mechanical forces/strains transmitting within the cytoskeleton. Bymodification of their configurations, the binding proteins activate upstreameffectors for the Rho GTPases. Once the Rho GTPases have been activated,the cell’s cytoskeleton starts to respond and produce effects such as motility.1.4 Experimental techniquesSo far, a wide variety of experimental devices have been applied to probethe cell’s mechanical properties [114, 205, 206]. Here, I only briefly reviewthe techniques closely related to this thesis.1.4.1 Micropipette aspirationIn micropipette aspiration experiments, a cell is aspirated into a mi-cropipette by applying a suction pressure. With a specific sucking pressure,the length of the tongue sucked into the pipette can be measured. By com-bining the experimental data with theoretical analyses [240] or numericalsimulations [5, 36, 73], one may determine the elastic and viscoelastic prop-erties of the cell. The micropipette aspiration method has been applied toinvestigate the heathy red blood cells [49], malaria infected red blood cells[62, 153], and white blood cells [51, 244], among other cell types.1.4.2 Microfluidic channelsA typical microfluidic channel has a narrow section connecting two reser-voirs. A pressure gradient is used to push the cells through the channel. Oneadvantage of using the microfluidic channel is that it mimics the in vivo en-vironment. This method can be used to measure the deformability of thecell, or distinguish and separate abnormal cells from healthy ones based ontheir different stiffness. The microfluidic channels have been used to testthe deformability of healthy red blood cells [66], malaria-infected red bloodcells [13, 67, 195], and neutrophils [237, 238].41.5. Cell models1.4.3 Substrate stretcherIn vivo, cells are often under cyclic stretch and compression. Examplesinclude vascular endothelial cells, airway smooth muscle cells and heart mus-cle cells. The substrate stretcher is usually used to study the response ofadherent cells, particularly the dynamics of the SFs and focal adhesions, un-der repeated stretch-compress cycles. Typically, the tested cells are culturedon an elastic substrate or membrane. Then, external mechanical forces areapplied to stretch or compress the substrate or membrane to deform theadherent cells. Such devices have been used to study cells under a singlestretch-compress maneuver [23, 59, 109, 213] and repeated stretch-compressmaneuvers [99].1.5 Cell modelsIn earlier cell modeling studies, the cell was treated as a continuum,either an elastic body or viscous droplet. The choice between these twodepends on the mechanical stiffness of the cell. Typically, suspended cellsbehave more or less as liquid drops, whereas adherent cells exhibit a higherelastic modulus and behave more like elastic bodies. These purely me-chanical models have contributed to measuring the cell’s viscous or elasticproperties in combination with experimental tools [29, 87]. Because of themechanosensing mechanism of the cells, chemo-mechanical models are re-quired to capture their dynamic response to mechanical perturbations. Inthis section, I first review the major mechanical models that have been pro-posed to study the cellular dynamics, including both continuum and discretemodels. Then, I will review some bio-chemo-mechanical models proposedfor studying the cell’s mechanosensing responses.1.5.1 Mechanical modelsThe purely mechanical models can be divided into two categories [118]:continuum models and discrete models.Continuum modelsSuspended cells have properties similar to droplets, examples being ery-throcytes and leukocytes in circulation. Thus, liquid-drop models have beenextensively used to study the rheology of suspended cells, e.g. neutrophilsentering into micropipette. The simplest fluid model is the Newtonian51.5. Cell modelsmodel, which treats the cell as a Newtonian droplet with a cytoplasmicviscosity. Such a simple model allows one to derivate analytical solutionsfor the cell dynamics under the micropipette aspiration. Yeung and Evans[240] and Needham and Hochmuth [154] derived constitutive equations forthe neutrophil flowing into the micropipette under an external sucking pres-sure. The model predicts a linear relationship between the sucking pressureand the rate of cell deformation, which agrees with some experimental ob-servations. Later, Tsai et al. [214] applied the Newtonian model to comparewith the micropipette aspiration experiments. They found that with highersucking pressure or flow rate they require a smaller apparent viscosity for thetheoretical model [154, 240] to predict the experimental measurements. Thisimplies a shear-thinning behavior of the neutrophil. Consequently, they pro-posed a power-law shear-thinning constitutive equation for the neutrophil.Drury and Dembo [42, 43] implemented the power-law shear-thinning consti-tution in their neutrophil model based on the finite element method. Theirnumerical simulations captured the acceleration of the cell at the end ofaspiration, in agreement with the experiment [214]. Later, Dong et al. [40]used a Maxwell constitutive equation to predict the elastic properties ofthe neutrophils at the beginning of the micropipette aspiration simulations.Despite these achievements, single-phase droplet models are too simple torepresent the complex inner structure of real cells. Later, Leong et al. [115],Marella and Udaykumar [128], and Zhou et al. [245] proposed compounddroplet models to study the passage of neutrophils in microfluidic channels.On the other hand, solid models are often applied to simulate cells withmuch higher stiffness, such as adherent cells [118]. The simplest one is thelinearly elastic model, which treats the cell as a homogeneous linearly elasticmaterial. Theret et al. [208] employed such a model to derive theoretical so-lutions for a cell under the micropipette aspiration. Later, Mijailovich et al.[136] used the linearly elastic model to investigate the responses of the cellin magnetic twisting cytometry experiments by finite element simulations.Viscoelastic models are also developed to capture the time-dependent me-chanics of the cell under complex load history. One intriguing example isthe linearly viscoelastic solid model that Sato et al. [184] used to analyzethe endothelium cell in the micropipette aspiration experiment.More complex mechanical models treat the cell plasma membrane as anelastic shell and the inner cytosol as a viscous liquid. I will call such modelselastic membrane models. Many numerical methods have been developedto simulate the elastic membrane models in simple flows. The plasma mem-brane is usually treated as a non-linear elastic shell. One widely used hyper-elastic constitutive equation for the membrane is Skalak’s model [199]. It61.5. Cell modelshas been used in simulation of capsules [6, 174], vesicles [239] and red bloodcells [90, 168, 204] under various flow conditions. An alternative consti-tutive equation for the hyperelastic membrane is the neo-Hookean model[6, 111, 112]. For the vesicle and red blood cell simulations, the bendingforce is accounted for. The most widely used bending model is the oneproposed by Helfrich [72, 190].Discrete modelsInstead of treating the cell components, such as the plasma membraneand the cytoskeleton, as continuum liquids or solids, discrete models em-ploy discrete particles to describe the cell geometry. Typically elastic springunits are used to link these particles to endow the cell with elastic resis-tance against deformation. Whereas continuum models rely on typicallyhomogeneous global constitutive equations to represent the cell’s mechani-cal properties, discrete models enjoy the advantage that the local materialproperties are represented by the local properties of particles and springs.Discher et al. [36] and Boey et al. [11] proposed an erythrocyte modelbased on a discrete particle-spring construct. In their model, each edgerepresents a spectrin polymer of the cytoskeleton of the red blood cell.This model has been applied to simulate the micropipette aspiration ex-periment [36]. One factor that hampers the application of such discretemodels in complex flows is the computational cost. Through coarse-grainingtechniques, the degrees of freedom of the discrete model can been signifi-cantly reduced. Such coarse-grained models include the dissipative parti-cle dynamics method [53–55, 165], the smoothed particle hydrodynamics(SPH) method [84, 85, 87] and the multiparticle collision dynamics method[155, 156]. These have found much application in studying the hydrody-namics of healthy and malaria-infected RBCs.Kamm and his coworkers et al. [12, 101, 102] have proposed a discretemodel to investigate mechanical properties of the cytoskeletal nework basedon Brownian dynamics. In their model, the actin filaments, binding pro-teins and myosin motors are modelled as discrete components separately.Their results indicated that all these components play crucial role in deter-mining the cytoskeletal properties. Such model may be further refined forinvestigating various actin-related phenomena.71.5. Cell models1.5.2 Chemo-mechanical modelsSome cells are able to respond to external mechanical forces by changingtheir cytoskeletal structures, which is termed cellular mechanosensing. Thecellular mechanosensing involves elements of external stimuli, biochemicalsignals, and mechanical response of the cells. So far, a few chemo-mechanicalmodels have been proposed to study the cellular mechanosensing.Models for adherent cellsInteraction between focal adhesions and stress fibers is a typical mechanosens-ing process relevant to most adherent cells. Novak et al. [157] proposed amathematical model to investigate the formation and distribution of the fo-cal adhesions under regulation of the contractile force. In the model, thecontractile force on the focal adhesions due to stress fibers is treated as aforce field. This model provides a mechanism to explain why more focaladhesions are located at peripheral regions of the adherent cells. However,it does not account for the stress fiber (SF) dynamics, which contribute tothe focal adhesion development.Models have been proposed to investigate SF dynamics under cyclicstretches. The SF is typically represented by a sarcomeric element andassumed to be tightly anchored to the substrate through focal adhesions.The sarcomeric model consists of an elasto-viscous element that connects toa contractile apparatus in series or in parallel. Kaunas and Deguchi [96] andHsu et al. [89] applied the sarcomeric model to simulate the SF dynamicsunder cyclic stretches, and investigated the effect of the myosin crossbridgecycling in regulating the SF tension. Later, Kaunas and Hsu [97], Kaunaset al. [98] and Tondon et al. [212] further applied the sarcomeric model tosimulate the reorientation dynamics of the SFs under cyclic stretches. Intheir simulations, assembly and disassembly dynamics of the SFs are mod-elled through a statistic kinetic equation. Stachowiak and O’Shaughnessy[202] used periodic sarcomeric model to study the SF retraction cut by lasernano-dissection. In all of these studies, the regulation of focal adhesionson the stress fiber is neglected. Chen et al. [17] performed an analyticalinvestigation of a sarcomeric model under cyclic stretch. By assuming acatch-bonding connection between the SFs and the focal adhesions, theypostulated that oscillatory forces shorten the life time of catch-bonds, whichcauses sliding of the focal adhesions and results in reorientation of the SFsperpendicular to the stretch direction. Later, Qian et al. [172] extended thistheoretical model into a 2D chemomechanical model that accounts for the81.5. Cell modelsfocal adhesion and SF dynamics, and gave a potential connection to the roleof Rho in the SF formation. They found that the reorientation dynamicsof the SFs under cyclic stretches depends on the competition between twofactors: the elevated forces within the SFs promoting disassembly and stiff-ening of the substrate stabilizing the focal adhesions. Besser and Schwarz[8] applied the sarcomeric model to simulate SF contractile dynamics underregulations of biochemical signal pathway. A positive feedback between theSFs and the focal adhesions is examined. The RhoA-pathway is simulatedby a reaction-diffusion equation based on the 1D stress fiber. This modelpredicts results that agree very well with experimental data.Using the continuum framework, alternative models have been proposedto handle the adherent cell mechanosensing [33–35, 225]. For example, thechemo-mechanical model of Deshpande et al. [33] assumes (1) a biochemicalsignal that regulates actin polymerization and myosin activation, (2) assem-bly of the actin and myosin into the SFs depending on the tension, and (3)contractility of the cell due to the SF tension through a cross-bridge cycling.The model successfully captures the experimental observations, such as highconcentrations of the SFs near the focal adhesions. Later, this model hasbeen extended to account for the SF disassembly [34] and the focal adhesiondynamics [35]. A similar continuum-based model is the force-dipole model,in which a force-dipole is used to represent a cell. De et al. [30] applied sucha model to predict the phenomena of cell reorientation under external cyclicstretch. The main drawback of this model is its difficulty in explaining thebiophysical mechanism for the cell reorientation.Models for suspension cellsHerant et al. [73–75] proposed a neutrophil model that accounts for dy-namics of the cytoskeleton that leads to a swelling force and a polymerizationforce. Using this model, they simulated the micropipette aspiration experi-ment. This work is particularly interesting to us because of the prominentroles that fluid dynamics plays in the process.9Chapter 2MethodologyThe three research projects to be presented rely on some common meth-ods for modelling and simulation. These include the smoothed particle hy-drodynamics (SPH) model for simulating fluid flow and solid-fluid interac-tion, the immersed boundary method (IBM) for handling no-slip boundaryconditions on complex fluid-solid interfaces, and various cell models. Theseare summarized in this chapter. The ensuing chapters, devoted to each ofthe three projects, will offer more specific details.2.1 Smoothed Particle Hydrodynamics (SPH)As a tool for computing fluid flow, the SPH method has been used widely,and several comprehensive reviews have appeared in recent years [144, 145].For our current purpose, we will only list the governing equations and out-line the solution algorithm. In the SPH formalism, continuum governingequations are discretized onto discrete particles more or less uniformly dis-tributed in space.2.1.1 Kernal integral representationThe concept of integral representation of a function A(r) in the SPHmethod starts fromA(r) =∫ΩA(r′)δ(r − r′)dr′ (2.1)where A is a function of the position vector r, and δ(r−r′) is the Dirac deltafunction. Ω represents the integral domain. The integral representation inequation (2.1) is always valid, as long as the A(r) is continuous in Ω.If the Delta function kernal δ(r−r′) is replaced by a smoothing functionW (r − r′, h), the integral representation of A(r) is rewritten as< A(r) >=∫hA(r′)W (r − r′, h)dr′ (2.2)102.1. Smoothed Particle Hydrodynamics (SPH)Figure 2.1: The support domain of a particle at r.where < · > represents the kernel approximation. W is the smoothingfunction. Figure 2.1 shows the diagram for the smoothing function at r.The smoothing function is usually chosen according to certain principles,which will be discussed in section 2.1.3. The error due to applying the SPHintegral representation can be estimated using the Taylor series expansionof A(r′) around r. Liu and Liu [120] proved that the kernal approximationof a function is of second order accuracy.The kernel approximation for the spatial derivative ∇·A(r) is obtainedby substituting A(r) with ∇·A(r) in the equation (2.2), which reads< ∇·A(r) > =∫h[∇·A(r′)]W (r − r′, h)dr′=∫h∇·[A(r′)W (r − r′, h)]dr′ −∫hA(r′) · ∇W (r − r′, h)dr′The first integral on the right hand side of the above equation can be trans-fered to integral over the surface S of the support domain. The smoothingfunction has zero values at the boundary of the support domain S, thus thefirst integral equals zero. We can rewriten the equation as< ∇·A(r) >= −∫hA(r′) · ∇W (r − r′, h)dr′. (2.3)112.1. Smoothed Particle Hydrodynamics (SPH)2.1.2 Particle approximationIn the SPH method, the whole computational domain is discretized by afinite number of particles that carry individual mass and occupy individualspace. If the element volume dr′ at the particle j is replaced by δVj =mjρj,then the continuous SPH integral representation for A(r) can be written inthe following form of discrete particle approximation< A(r) >=∫hA(r′)W (r − r′, h)dr′ ≈∑hmjρjA(rj)W (r − rj , h) (2.4)where∑h denotes the summation of particle at the support domain h.The particle approximation of function A and gradient of function A atthe particle i can be written as< A(ri) > ≈∑hmjρjA(rj)Wij , (2.5)< ∇·A(ri) > ≈∑hmjρjA(rj) · ∇iWij (2.6)where Wij = W (ri − rj , h). The gradient ∇i ·Wij is taken with respect tothe particle i:∇iWij =ri − rj|rij |∂Wij∂rij(2.7)There are a number of ways to derive SPH formulation of PDEs. In orderto obtain a symmetric form of particle approximation for gradient operators,Monaghan [142] applied the following form∇·A(r) =1ρ[∇·(ρA(r)−A(r) · ∇ρ] (2.8)∇·A(r) = ρ[∇·(A(r)ρ2) +A(r)ρ2· ∇ρ]. (2.9)The particle approximation for the above equations can be written as< ∇·A(ri) >≈1ρi[∑hmj [A(rj)−A(ri)] · ∇iWij ] (2.10)< ∇·A(ri) >≈ ρi[∑hmj [A(rj)ρ2j+A(ri)ρ2i] · ∇iWij ]. (2.11)122.1. Smoothed Particle Hydrodynamics (SPH)2.1.3 Kernel smoothing functionLiu and Liu [120] disccused a list of conditions that one should obey inconstructing a kernel smoothing function:• The normalization condition∫hW (r − r′, h)dr′ = 1where h represents the support domain of the smoothing function;• The Delta function limitation when the smoothing length approcheszerolimh→0W (r − r′, h) = δ(r − r′)• The compact conditionW (r − r′, h) = 0, when|r − r′| > κhwhere κ is a constant related to the smoothing function for point at r,and defines the effective support domain of the smoothing function.• The positivity conditionW (r − r′, h) > 0,for any point at r′ within the supporting domain.• The decay conditionThe function monotonically decreases with the increase of r′.• The symmetry conditionW (r − r′, h) = W (r′ − r, h),for any point at r′ within the supporting domain.• The smoothness condition: the kernal function and its derivative shouldbe smooth in the whole domain.So far, many forms of kernel smoothing functions have been proposed[120]. In this section, I will only display the quadratic smoothing function132.1. Smoothed Particle Hydrodynamics (SPH)proposed by Johnson et al. [93], which has been applied in my RBC simu-lation. The quadratic smoothing function readsW (r, h) = αd(316(rh)2 −34rh+34), 0 ≤ R ≤ 2h (2.12)where αd is the dimension-related coefficient, which equal to 2pih2 for 2D and54pih3 for 3D respectively.2.1.4 Weakly compressible SPH algorithmIt is possible to apply projection schemes to impose incompressible con-dition in the SPH method [86, 113, 193]. But we have adopted a weaklycompressible scheme in the current research for simplicity. The detailedscheme for advancing the particle motion will be discussed in this section.We have the following equations governing the motion of fluid particles:DvDt= −1ρ∇P +1ρ∇ · τ + fm, (2.13)DρDt= −ρ∇ · v, (2.14)PP0=(ρρ0)7− 1, (2.15)where D/Dt is the material derivative taken on each moving particle, v isthe velocity of the particle, ρ is the density of the particle, P and P0 are theinstantaneous pressure and the reference pressure of the particle. τ is theviscous stress tensor. fm is the membrane force acting on the membraneparticles. For fluid particles, this force equals zero. The details of how toderive this term will be discussed later.The momentum equationThe pressure gradient term was used a symmetric form as∇Piρi= −∑hmj(Piρ2i+Pjρ2j)∇iWij .The laminar stress formula [194] is applied, which reads1ρi∇·τ i =∑h4ν0mjrij(ρi + ρj)|rij |2vij∇iWij ,142.1. Smoothed Particle Hydrodynamics (SPH)where ν0 is the kinetic viscosity of the fluid. Putting the above equationsinto the momentum equation (2.13), we can getDviDt=−∑hmj(Piρ2i+Pjρ2j)∇iWij+∑hmj4ν0rij(ρi + ρj)|rij |2vij∇iWij + f im. (2.16)The continuity equationIn the SPH, one extensively applied method to estimate the density ofthe particle i is the summation density approach [142], which readsρi =∑hmjWij (2.17)An alternative is to track the evolution of the particle density through amass conservation equation, which can be written in particle representationasDρiDt=∑hmj(vi − vj)∇iWij . (2.18)Two advantages can be gained by applying the latter form (equation2.18). (1) The density variation only depends on relative motion betweenparticles, which avoids density oscilation on the particles near the edge byusing the equation (2.17). (2) The rates of change of all physical variablescan be computed in one subroutine, which saves computation as comparedwith doing summation through the equation (2.17).The equation of stateIn compressible flows, the density variation δρ/ρ0 is proportional to M2,where M is the Mach number for the local flow. If the fluctuation of thedensity is required to be smaller than η in one problem for instance, theMach number should satisfy the condition ofM =vfcs<√ηwhere vf is the local fluid velocity. cs is the sound speed of the flow, whichcan be calculated by cs =√∂P/∂ρ|ρ0 . In the current study, the pressureand the density are related through the Tait equation [143]P = P0[(ρρ0)γ − 1] (2.19)152.2. Immersed Boundary Methodwhere γ is typically set to be 7. Based on the above equation, we canobtain the formula for the local sound speed: cs =√P0γρ0. The condition forconstraining the density variation less than η isP0 >ρ0vfγη. (2.20)where vf is the local fluid velocity. In my simulation, P0 is choosen that thedensity deviation for all fluid particles is less than 1%.2.1.5 Boundary treatmentIn the SPH method, the particles near the boundary suffer from the issueof truncated support domains [120]. The implementation of the boundaryconditions for the SPH method is not as simple as in the mesh-based meth-ods.For the wall boundary condition, several techniques have been proposedto solve the insufficiency of the particle support domain for the particlesnear the wall.• ghost particles: Monaghan [143] proposed to apply a line of “ghost”particles to represent the wall boundary, which produce repulsive forceson the approching particles to avoid penatration of the particles throughthe wall.• image particles: Libersky et al. [117] introduced multi-layer imaginaryparticles to reflect a symmetrical surface boundary condition with op-posite velocity on the reflecting image particles.In the iRBC simulation of this thesis, “image particles” are employed toimpose the non-slipping wall boundary condition. In the flow direction, Iapplied the periodic boundary condition, which assigns the particles thatexist the channel back to the entry end of the channel.2.2 Immersed Boundary Method2.2.1 Governing equationsThe fluid dynamics is governed by imcompressible Navier-Stokes equa-tions:∂v∂t+ v · ∇v = −1ρ∇p+µρ∇2v + fw + fm (2.21)∇·v = 0 (2.22)162.2. Immersed Boundary MethodFigure 2.2: (a) Diagram for a solid body immersed in a flow. The bodytakes the volume Ωb with boundary Γb. (b) Schematic of body immersed ina Cartesian grid. From Mittal and Laccarino [138] with permission, c©2005,Annual Reviews.where v is the velocity, ρ is the density, and µ is the viscosity of the fluid.fw and fm are force acting on the Eulerian nodes due to immersed solidwall and membrane, respectively.2.2.2 Imposition of immersed boundariesWe treat the fluid-solid interaction by the immersed boundary method[52, 163]. In the current context, we use capital letters (i.e., Fw and Fm) todenote quantities on the Lagrangian nodes, and lowercase letters for thoseon the Eulerian nodes (i.e., fw and fm). Transformation between the La-grangian and Eulerian quantities is realized through interpolation.Elastic membrane boundariesThe treatment of the elastic membrane follows that of Peskin [163]. Theelastic membrane is represented by a collection of massless points connectedby springs. The coordinate Xi of the ith Lagrangian point is governed bythe equation∂Xi∂t= v(Xi, t), (2.23)where v(Xi, t) is the velocity of the ith Lagrangian point. Deformation ofthe membrane configuration (i.e., Xi) produces a membrane tension. This172.2. Immersed Boundary Methodforce can be distributed to the surrounding Eulerian nodes byf(x, t) =∑iF i(t)δ(|x−Xi|), (2.24)where δ is the smoothed Dirac delta function, and F i(t) is the elastic forceacting on the ith Lagrangian node. In this study, it refers to the force onthe ith membrane particle, i.e., F im in the Eq. (4.6). The Lagrangian pointsdo not necessarily coincide with the Cartesian nodes, of course. Thus asmoothed distribution function is used to replace the sharp delta function.Rigid wall boundariesOur simulation will also involve solid boundaries in contact with fluid,at which the no-slip condition must be enforced. We adopt the so-calleddirect forcing method for this purpose [52]. This can be illustrated by thediscretized momentum equation for the fluid flow:vn+1 − vn∆t= RHS + fw (2.25)where RHS is the discretized form of [−v · ∇v − ∇pρ +µρ∇2v]. Its actualform depends on the computational scheme used, and is immaterial to ourpurpose here. The “wall force” fw vanishes on Eulerian grid points insidethe fluid bulk, and will be chosen for wall grids so as to enforce the no-slipboundary condition: v = V w. By replacing vn+1 by V w, we getfw =V w − vn∆t−RHS. (2.26)This determines the “direct forcing” fw that is required in Eq. (2.25) sothat the fluid velocity equals the wall velocity V w.2.2.3 The smoothed δh functionSince the Lagrangian nodes are usually not coinsident with the Euleriannodes, a smoothed δh function is required to distribute the Lagrangian valuesto the Eulerian nodes, as is shown in the Fig. 2.3. The three dimensionalsmoothed δh function is usually written asδh(r) =1h3φ(rxh)φ(ryh)φ(rzh) (2.27)182.2. Immersed Boundary MethodFigure 2.3: (a) Distribution of forcing F k from Lagrangian boundary point(Xk) to surrounding fluid nodes. Shaded region denotes the force distribu-tion region. (b) Distribution functions employed in various studies. FromMittal and Laccarino [138] with permission, c©2005, Annual Reviews.where h is the resolution of the Eulerian mesh, and rx, ry and rz are Carte-sian components of r. We employed the φ function proposed by Peskin [163]in the current study, which readsφ(x) =0, x ≤ - 2(5 + 2x−√−7− 12x− 4x2)/8, -2 ≤ x ≤ -1(3 + 2x+√1− 4x− 4x2)/8, -1 ≤ x ≤ 0(3− 2x+√1 + 4x− 4x2)/8, 0 ≤ x ≤ 1(5− 2x−√−7 + 12x− 4x2)/8, 1 ≤ x ≤ 20, 2 ≤ x2.2.4 Computational algorithmThe fluid subproblem is a standard Navier-Stokes problem. The readeris referred to Yu et al. [241, 242] for the detailed description of our solver.As an example, the momentum equation is discretized asv∗ − vn∆t−µ2ρ∇2v∗ =−∇p−12[3(v · ∇v)n − (v · ∇v)n−1]+µ2ρ∇2vn + fnw + fnm, (2.28)∇·v = 0. (2.29)192.3. Cell models2.3 Cell models2.3.1 Red blood cell (RBC) modelThe cell membrane is represented by a discrete triangular network, whichis endowed with elasticity against in-plane strain, as well as bending. Be-sides, energetic constraints are applied to maintain the cell volume and sur-face area. The cell is immersed in the discrete fluid particles with the mem-brane particles marked. The membrane forces are applied directly on themarked membrane particles. The detailed forms for these forces are givenin section Neutrophil cell modelThe neutrophil membrane is represented by a discrete spring-particlecircle, which is endowed with elasticity against in-plane strain, as well asbending. Besides, energetic constraints are applied to maintain the cell vol-ume and surface area. The cell is immersed in the continuum fluid medium.The elastic membrane forces are distributed to surrounding fluid mediumby the immersed boundary technique. The detailed forms for these forcesare given in section Stress fiber model for adherent cellsIn the substrate stretcher experiment, the mechanical properties of theadherent cells are reflected by the stress fibers (SFs) that are anchored to thesubstrate through focal adhesions. In this case, the fluid hydrodynmics donot play significant role in the cellular dynamics. Besides, Trepat et al. [213]and Krishnan et al. [109] have revealed that fluidization of the cell is due todisassembly of the SFs. Thus, I only focus on studying the dynamics of singleSF. In the current study, the SF is represented by a mechanical Kelvin-Voigtelement connected to a myosin apparatus. Please refer to section 5.2 for thedetailed formulation of the SF model.20Chapter 3Simulation ofmalaria-infected RBC inmicrofluidic channels:Passage and blockage3.1 IntroductionHealthy red blood cells (RBCs) have a biconcave shape with a diame-ter of 7–8 µm and a thickness of 2–3 µm [50]. During their normal func-tion, RBCs traverse the circulatory system, squeezing through capillariesas small as 3 µm in diameter. Their high deformability is mainly due totheir excess surface area and extremely flexible membrane. When infectedby the malaria-causing parasite Plasmodium falciparum, however, the RBCgradually loses its deformability through the three stages of infection—ring,trophozoite and schizont—as the parasite grows and multiplies asexuallywithin the infected RBC (iRBC) [27, 71, 141]. A direct consequence of thehardening of the iRBC is occlusion of microcapillaries, which causes themost severe symptoms of malaria [137].The rigidification of iRBC has been quantified using several experimen-tal protocols: micropipette aspiration [153], stretching by optical tweezers[29, 55, 87], and passage through microfluidic channels [67, 76, 77, 195].Aspiration and stretching of RBCs directly test their deformation undermechanical forcing, and have yielded considerable insight into the harden-ing of the cells. So far, three mechanisms have been identified for the lossof deformability of the iRBC. First, the iRBC membrane hardens as par-asitogenic proteins induce abnormal cross-linking of the spectrin network[62, 63]. Second, transport of liquid and proteins across the compromisedcell membrane alters the iRBC’s surface-to-volume ratio and swells the cell[47, 192]. Finally, the solid parasites inside the iRBC, known as merozoites,mechanically hinder cell deformation, especially when they grow in size and213.1. IntroductionFigure 3.1: Passage and blockage of microfluidic channels of different sizesby iRBC in the ring, trophozoite and schizont stages of infection. The flowgoes from the right to the left. From Shelby et al. [195] with permission,c©2003 The National Academy of Sciences.multiply by mitosis in the schizont stage [87, 153].With the advent of microfluidic assaying, microfluidic channels have beenincreasingly used as a tool for probing the mechanical property of cells[13, 88, 217]. For malaria-infected RBCs, microfluidics brings one greatadvantage: it mimics the geometry in microcapillaries and in the spleen[16], where the red cells are forced through narrow passages. Consequently,it may potentially bridge the gap between idealized mechanical measure-ments using micropipettes and optical tweezers and in vivo events of interest,such as occlusion of capillaries and filtering of hardened RBCs in the spleen[16, 27, 38, 186]. Several groups [67, 76, 77, 195] have employed microflu-idic channels to probe the deformability of healthy and malaria-infectedred cells. Typically, a prescribed pressure drop is applied to push iRBCsthrough microfluidic channels of various sizes. For example, Shelby et al.[195] demonstrated that with progression of the infection, the iRBC hard-ens gradually, and blocks larger and larger pores (Fig. 3.1). So far, however,experimental work in this area has remained largely qualitative. It would bedesirable to establish a quantitative connection between the geometric andmechanical attributes of iRBCs and their ability to traverse microchannels223.2. Physical model and numerical schemeor blood vessels of a certain size. In the development of biomechanical assay,Ma and coworkers [67, 67, 134] measured the critical pressure required topush a cell through a contraction, and extracted a cortical tension throughthe Young-Laplace equation. This only serves as a general indicator of thecell’s rigidity. However, it does not account for the non-uniform deforma-tion of the cell membrane, and thus cannot be quantitatively related to themechanical and structural changes in the iRBC.This study aims at establishing a quantitative relationship between cellrigidity and passage or blockage of microchannels. By using a particle-based formalism in three dimensions (3D), we simulate the transit of healthyand infected RBCs through a prototypical microfluidic channel with a con-traction. The surrounding fluid and the red cell, including the membrane,cytosol, and possibly the merozoite, are all discretized by particles in theframework of smoothed particle hydrodynamics [87]. The interaction amongthe membrane particles is designed to produce an overall constitutive be-havior similar to the commonly used hyperelastic models. By studying thetransit time and critical condition for occlusion, we delineate the contribu-tion to iRBC rigidification from each of the three mechanisms, hardening ofthe membrane, loss of excess surface area, and growth of the merozoite. Theresults not only shed light on in vivo RBC behavior in microcapillaries andin the spleen, but also provides guidelines for designing sensitive microfluidicassays.3.2 Physical model and numerical schemeMost of the cell-mechanics simulations so far are based on continuummodels, ranging from the liquid drop model to vesicle models [115, 118, 187–189, 240]. Particle-based models are relatively new, and by virtue of theirmeshless formalism, enjoy greater convenience in representing morphologicaland structural changes in the cell or intracellular components [28, 53–55, 87,91, 160, 165, 215]. It is mostly this advantage that drew us to a particle-based model.The physical model in this work is similar to that used by Hosseini andFeng [87]. We discretize the cytosol and surrounding plasma by particlesmuch as in traditional Smoothed Particle Hydrodynamics (SPH) modelsfor flow computation. The cell membrane is tesellated into more or lessregular triangles, the vertices being represented by particles and the edges bysprings. Thus, the membrane is discretized into a particle-spring network.The parasite is realized numerically by constraining a set of particles to233.2. Physical model and numerical schememove as a rigid body [144], which floats in the cytosol. The motion of the“fluid particles”, representing the cytosol and plasma, is computed usingthe conventional SPH method and reproduces the fluid flow governed bythe Navier-Stokes equations. The membrane particles experience bendingand in-plane elasticity in additional to fluid pressure and viscous stress, asexplained below.The numerical solution is based on the open-source SPH solver parallelSPHysics [176], developed for simulating fluid flows with free surfaces [64,65]. On the basis of the original parallelized FORTRAN code, we have addedthe treatment of the cell membrane. The parallelization is realized throughdomain decomposition, and the interaction among particles that reside indifferent cores are realized by the MPI formalism [176].3.2.1 Cell membrane elasticityThe cell membrane is represented by a discrete particle-spring network,which will be endowed with elasticity against in-plane strain and bending.We use linear springs with a stretching coefficient of ks such that the elasticenergy for in-plane deformation isEs =∑i,jks2(Lij − Lij,0)2, (3.1)where the summation is over all pairs of adjacent vertices i and j, Lij isthe length of the spring connecting them, and Lij,0 is its resting length. Forbending, Hosseini and Feng [87] applied the Helfrich bending energy [190]to compute the distribution of bending forces on the membrane particles.We found this formalism to be quite sensitive to the quality of the mesh;when the cell undergoes large deformation, the distorted triangles tend toproduce elastic instability due to inaccurate evaluation of the local mem-brane curvature. Thus, we decided to adopt the simpler and more robustbending model of Wada et al. [216, 224], with a bending coefficient kb andbending energyEb =∑i,j2kb tan2(θij2), (3.2)where the summation is over all pairs of neighboring triangles i and j and θijis the angel between their normals. Note that this bending energy assumeszero spontaneous curvature for the membrane.In addition, the red cells are known to conserve their surface area, and inour particle model, this is implemented through an energy penalty against243.2. Physical model and numerical schemelocal area dilatation:EA =kd2N∑j(Aj0 −AjAj0)2Aj0, (3.3)where kd is a constant, Aj0 is the undeformed area of the jth triangle, andthe summation is over all N triangles of the RBC membrane. Finally, weinclude an energy penalty against the change of the total cell volume:EV =kv2V0(V0 − VV0)2, (3.4)where kv is a constant coefficient, and V0 is the initial volume of the cell. Un-der large forcing and severe cell deformation, this volume constraint helps toprevent fluid particles from penetrating the membrane. The two constrainsin Eqs. (3.3) and (3.4) are similar to those used by Discher et al. [36] andFedosov et al. [55], but we have excluded their constraint on the total sur-face area. In all simulations, the cell surface area and volume are conservedto within 2%. Using Eqs. (3.1–3.4), we write the total elastic energy of thecell membrane as Em = Es +Eb +EA +EV . The elastic force acting on themembrane particles can then be calculated asfm = −∂Em/∂r. (3.5)3.2.2 Parameters for the modelThe elasticity of our discrete network of particles and springs can berelated to continuum models for an elastic membrane by considering thelimit of small strain. For a regular and isotropic triangular mesh, the net-work has an effective shear modulus Gs =√34 ks, Young’s modulus E =2√3ks, Poisson ratio ν = 13 and area dilation modulus K =√32 ks [160].The bending constant kb can be related to the Helfrich coefficients as [9]:kb = 2√3(κc + κg/2) = 2√3κ, κc and κg being the coefficients for the meancurvature and the Gaussian curvature terms in the Helfrich bending energy,and κ the average bending modulus [36, 55]. There is a small amount ofnonuniformity in our tessellation, necessitated by the requirement of cov-ering the curved biconcave surface of the RBC. This may have introducedsmall deviations from the above relationships. For nonlinear elasticity underlarge strain, Omori et al. [160] have compared the strain-hardening of springnetworks with that of the continuum Skalak and Mooney-Rivlin models. For253.2. Physical model and numerical schemeTable 3.1: Dimensions of the healthy and infected red cell in different = S/(4pia2) denotes the excess surface area ratio of the cell, where a isthe effective cell radius defined as the radius of a sphere having the volumeV of the cell. The last column indicates the dimensions of the parasite.Case S (µm2) V (µm3) se a (µm) Parasite (µm)Healthy 132 92 1.34 2.80 NARing 132 92 1.34 2.80 R = 1, h = 0.5Early trophozoite 132 105 1.23 2.92 R = 1.4Late trophozoite 132 116 1.15 3.02 R = 1.7Schizont 97 85 1.04 2.72 R = 2.4, h = 3most of the simulations reported here, we have used a fixed bending constantkb = 6.93× 10−19 J, chosen according to Hochmuth et al. [81]. To examinethe role of in-plane elasticity, we will vary the stretching constant ks in arange consistent with experimentally measured membrane shear modulusin healthy RBC and iRBC in different stages of malaria infection. Basedon such parameters, the ratio between the bending forces and the in-planestretching forces can be estimated as ξ = kb/(a2ks) = O(10−3), where a isthe effective cell radius. Thus, bending contributes little to the membranedynamics and cell deformation.3.2.3 Cell geometriesFor the healthy RBC, we adopt the biconcave shape specified by theformula of Evans and Fung [50]:D(r) =√1− (r/R0)2 [C0 + C1(r/R0)2 + C2(r/R0)4], (3.6)where D(r) is the thickness of the RBC as a function of distance r from thecenter axis, and R0 is the radius of the RBC, measured from the center tothe outer edge (Fig. 3.2). We have taken the following values: R0 = 3.9 µm,C0 = 0.81 µm, C1 = 7.83 µm, C2 = −4.39 µm; the center has a thickness ofT1 = 0.81 µm and the largest thickness at the rim is T2 = 2.4 µm.Through the ring, trophozoite and schizont stages, the iRBC changes itsvolume, shape and surface area, while the parasite grows in size. A sur-vey of experimental measurements of iRBC and parasite dimensions revealsconsiderable variations not only among different studies, but also amongcells of the same cohort [47, 48, 153, 192]. We have strived to capture theoverall trend of these data, and adopted the dimensions in Table 3.1. In the263.2. Physical model and numerical scheme-440-2 20 24-2-4012-1-2R0T1 T2Figure 3.2: The healthy RBC and the ring-stage iRBC have the same bicon-cave shape described by Eq. (3.6), with the smallest and largest thicknessbeing T1 = 0.81 µm and T2 = 2.4 µm. The spatial coordinates are markedin microns.ring-stage, the shape and size of the iRBC do not differ appreciably fromtheir healthy counterparts [192], and we have adopted the shape in Fig. 3.2for the ring-stage iRBC as well. The merozoite is a disc of radius R = 1 µmand thickness h = 0.5 µm. Drastic morphological changes of the iRBC oc-cur in the trophozoite and schizont stages. During the trophozoite stage,the merozoite continuously digests cytoplasmic hemoglobin and grows from2.0 µm to 4.0 µm in diameter. Meanwhile, the cell surface area remainswithin a narrow margin of variation [47, 48, 153, 192]. Thus, the growingparasite progressively swells the cell and gives it an asymmetric biconcaveshape. Here we represent the iRBCs at the early and late trophozoite stagesby the shapes of Fig. 3.3, with the merozoite being a sphere of increasingradius R. In the schizont stage, the iRBC undergoes a marked reductionin volume and surface area [47, 161, 192, 196], and assumes a more or lessspherical shape [196]. The merozoites multiply and form a tightly packedcluster that may occupy as much as 80% of the iRBC’s volume [71, 192].We represent the schizont-stage iRBC as a spheroid with radius of 3.2 µmand thickness of 4 µm (Fig. 3.4). The cluster of merozoites is modeled as asingle spheroidal particle of the same aspect ratio as the iRBC, with radiusR = 2.4 µm and thickness h = 3 µm.The bending energy introduces a subtlety in realizing the prescribed cell273.2. Physical model and numerical scheme(a)-440-2 20 24-2-4012-1-2D1D2T2 T1(b)T2 T1D2D1Figure 3.3: The cell geometries at (a) the early trophozoite stage, with(D1, D2, T1, T2) = (7.77, 7.76, 2.7, 3.0) µm, and (b) the late trophozoitestage, with (D1, D2, T1, T2) = (7.80, 7.72, 2.8, 4.0) µm.shapes in our particle-spring model. At the start of a simulation, we initializethe resting length of all springs to that corresponding to the tessellation ofthe underformed cell surface. This relieved all in-plane stretching. Forthe healthy and ring-stage cells, the shape of Eq. (3.6) is very close to aminimizer of the bending energy [10, 190]. Our simulation indicated that inthe static fluid medium turning on the bending energy slightly derives thecell from its initial configuration to a new equilibrium one. The difference isless than 2% in the cell dimensions. For the trophozoite and schizont stages,on the other hand, the chosen iRBC geometries differ considerably from283.3. Results-440-2 20 24-2-4012-1-2DTFigure 3.4: The iRBC has a spheroidal shape at the schizont stage, withdiameter D = 6.4 µm and thickness T = 4 µm.the bending-energy minimizer. This does not constitute a complication inpractice because cell deformation will be dominated by hydrodynamic forcesand in-plane elasticity in our simulations, and bending has negligible effects.3.2.4 Repulsive force between membrane and wallFinally, as is commonly done in the literature [39, 115], we have alsoemployed a repulsive potential to avoid solid-solid contact. Under large de-formation, the cell membrane particles tend to penetrate into the wall or theparasite, e.g. when the cell squeezes through narrow channels. We adopt thecut-off spring model [39, 115], which uses a repulsive potential with a finiterange and a linear elastic constant kr. In our case, numerical experimentsshow that setting the range to the initial nearest-neighbor separation is agood choice. We have also tested different kr values, and found that theresults are insensitive to kr within the range of 2× 10−4 to 10−2 N/m. Wehave used kr = 2× 10−3 N/m for all the results reported.3.3 ResultsThe computational domain is schematically depicted in Fig. 3.5. Thegeometry and dimensions of the microfluidic channel are inspired by recentexperimental setups [67, 134, 195]. The entrance has a rectangular cross293.3. ResultsLlllwWHFigure 3.5: Schematic of the computational domain. The overall dimensionsof the channel are L = 24 µm, W = 12 µm and H = 4.8 µm, with threesegments of equal length l = 8 µm. The narrow section in the middle has awidth w = 3.2, 4 and 4.8 µm for the narrow, medium and wide channels.section of width W = 12 µm and height H = 4.8 µm. The same height ismaintained throughout the entire conduit, but the width contracts througha 45◦ shoulder to a narrower width w. Further downstream is a suddenexpansion to a cross-section that is identical to the one at the entrance.The length of the narrow segment is fixed at l = 8 µm, but its widthtakes on one of three values, w = 3.2, 4 and 4.8 µm, which will be calledthe narrow, medium and wide channel, respectively. The red cell alwaysenters “edgewise”, that is, with its circular mid-plane horizontal and midwaybetween the top and the bottom of the microchannel in the schematic.Based on physiological values [79, 100], we have set the density and vis-cosity of the blood plasma (or surrounding fluid in the microfluidic channel)to be µ = 10−3 Pa·s and ρ = 103 kg/m−3. In healthy red cells, the cytosolviscosity is typically 3–6 folds higher than the plasma viscosity because ofthe hemoglobin and other proteins in the cytosol [20]. This viscosity ratiocan affect the cellular dynamics [25]. Moreover, during malaria development,the parasite continually digests the hemoglobin. Thus the cytosol viscositymay evolve in time, although no quantitative results have been found in the303.3. Resultsliterature. For simplicity, we have assigned equal viscosity to the cytosoland the plasma (λ = 1) as did many previous simulations [53, 204]. Besides,the density of all fluid and solid components is set to ρ = 103 kg/m−3.To simulate a pressure-driven flow through the microchannel, we imposeperiodic boundary conditions between the entrance and the exit of the chan-nel, and apply a constant body force on all the particles to drive the flow.The body force per unit mass Fx can be related to the effective pressure dropacross the length of the channel as such: ∆P = ρFxL. This setup is prefer-able to imposing a pressure gradient directly because the SPH methodologyhandles the periodic boundary conditions most easily [86]. From this pres-sure drop, we can define a characteristic velocity U = ∆Pw2/(8µL) and acharacteristic time T = L/U . Now a capillary number can be defined:Ca =µUGs. (3.7)In the results to be presented below, we fix ∆P to be 12 Pa, on the sameorder of magnitude as in recent microfluidic experiments on iRBC deforma-bility [13, 67]. Under these conditions, Ca lies between 2.56 × 10−3 and2.88 × 10−2 for the range of membrane modulus Gs to be examined. TheReynolds number, defined using the characteristic velocity U , is on the orderof 10−2.To establish numerical convergence of the results with respect to spatialresolution, we have simulated passage of the healthy RBC through the widechannel (w = 4.8 µm) using two resolutions. The coarse resolution has anominal particle separation l0 = 0.6 µm, with 440 particles on the RBCmembrane and 10977 particles overall. For the fine resolution, l0 = 0.4 µmand there are 1058 particles on the membrane and 29648 particles overall.The cell trajectories for these two cases agree within 3% throughout thepassage, and the overall transit time differs by less than 5%. We havefurther checked the cell shape and the membrane strain contours, and thedifferences between the two cases are very small. All subsequent results arepresented using the fine resolution.As explained in the Introduction, the rigidification of iRBCs has threepotential causes, the stiffness of the membrane, the cell geometry, and theinternal components, including the shape and size of the merozoite. In thefollowing, we will first investigate each element separately while fixing theother parameters at baseline values. Then we will combine all three changestogether as occurs during a real infection, and examine how the overall rigid-ity of the iRBC evolves through the ring, trophozoite and schizont stages.313.3. Results3.3.1 Effect of the membrane stiffnessAccording to prior measurements using micropipette aspiration, the healthyRBC has a membrane shear modulus around 5 µN/m [50]. After infectionby malaria, the merozoites export parasitogenic proteins onto the membranethat remodel the spectrin network [125, 181]. A direct consequence of thismodification is an increase in the membrane rigidity, with the shear modu-lus Gs being elevated by 3 to 5 folds in the trophozoite and schizont stages[63, 87, 153]. In this subsection, we fix the resting shape of the cell to that ofthe healthy RBC, and systematically examine the influence of increasing theshear modulus Gs, from 5 µN/m to 35 µN/m. No merozoite will be includedinside the cell. Specifically, we will study the passage or arrest of the cells inwide, medium and narrow microchannels. The change in Gs is effected byvarying the stretch constant ks of the springs in our particle-spring network.The Young’s modulus and area dilation modulus will change simultaneouslyas well.Figure 3.6 depicts the transit of the RBC through the narrow channel atthree values of the shear modulus. Similar trends are seen in the mediumand wide channels. The instantaneous position of the RBC is indicated byxc, the longitudinal coordinate of the centroid of all membrane particles,scaled by the channel length L. Time is scaled by T . The horizontal dashedlines indicate the entry and exit to the narrow middle segment. Since theundeformed RBC has a diameter of 7.8 µm, more than twice the width ofthe narrow passage, the cell undergoes large deformations. The cell with thesoftest membrane (Gs = 5 µN/m) passes through the contraction readily,while the one with the intermediate Gs takes longer time. The cell with thehighest membrane modulus (Gs = 25 µN/m) fails to traverse the channel; itblocks the entry to the narrow segment. The insets show snapshots for thelatter two runs. On the surface of the cell, we have superimposed contourplots of the local deformation γ. For the ith membrane particle, γi is definedby averaging the deformation of all the Ni springs connected to it:γi =√√√√∑Nij=1(LijLij,0− 1)2Ni, (3.8)where Lij and Lij,0 are the current and resting lengths of the springs, and Niis the number of springs linking to the ith membrane particle. Generally, themembrane suffers large shear deformation in the dimple region, i.e. the thincentral area of the undeformed RBC, and relatively small strains at the rim.The stretching γ can reach 30% for Gs = 15 µN/m after the cell enters the323.3. Results0. 3.6: The trajectories of RBCs through the narrow channel (w =3.2 µm) for three values of the shear modulus Gs. The insets show top-view snapshots of the cell during the transit with color contours of thelocal stretching γ of Eq. (3.8). The red arrow pointed insets are for Gs =15 µN/m at t = 0.27 and t = 1.9, while the blue arrow pointed inset is forGs = 25 µN/m at t = 2.1.narrow channel completely. Also, we notice that the squeezing has producedone or more longitudinal wrinkles on the membrane, typically at the centraldimple. Because of the slight asymmetry in the initial tessellation of thesurface, the wrinkle and the stretching contour appear slightly asymmetric.The most obvious effect of Gs is in slowing down the transit. This isevident from comparing the trajectories for Gs = 5 µN/m and 15 µN/m inFig. 3.6. To be more quantitatively, we define a transit time tt as the intervalthat starts with the cell’s forefront reaching the entry and ends with the itscentroid leaving the exit. Figure 3.7 plots the transit time as a functionof the membrane shear modulus for the three channels. Membrane rigidityis seen to lengthen the transit time in all channels. Moreover, this effectbecomes more pronounced for narrower channels. For the wide and mediumchannels (w = 4.8 and 4 µm), increasing Gs from 5 to 35 µN/m increases tt333.3. ResultsFigure 3.7: The transit time tt, made dimensionless by T , as a function ofthe membrane shear modulus Gs for the three about 20% and 30%, respectively. For the narrow channel (w = 3.2 µm),on the other hand, tt diverges at Gs = 25 µN/m; the cell blocks the channeland fails to pass.Cell blockage at Gs = 25 µN/m indicates that the supplied ∆P , cor-responding to Ca = 2.56 × 10−3 in this case, is insufficient to deform thecell to the extent required for passage. Numerical experiments show thatdoubling ∆P or Ca will lead to passage of this cell. Thus, there exists a crit-ical pressure for passage between these two situations, and such a pressureis expected to correlate directly with the membrane modulus. Previously,Ma and coworkers [67, 134] have used such a critical pressure to probe thecortical tension of red blood cells. This tension is assumed to be constantin the membrane, regardless of the local deformation, and is extracted fromthe critical pressure drop via the Young-Laplace equation. With the helpof simulations that properly account for the membrane elasticity, a similarprocedure can be developed as a means of measuring the cell membranemodulus.We end this subsection by pointing out a subtle feature of the blockage.In microfluidic experiments, blockage is often incomplete. Because of the343.3. Resultsbending elasticity of the membrane, the corners of the typically rectangularcross-section may not be blocked, and a fluid flow persists through thesesmall openings [61, 67, 195]. Our simulations appear also to show a thingap between the wall and the cell membrane; see insets of Fig. 3.6. Thisclearance is due to the repulsive force between the cell membrane and thewall, and is roughly equal to the nominal particle separation l0. Thus,it is too small for particles to slip through, even at the corners, and theblockage of the microchannel is complete. This subtle difference betweenexperimental and computational blockage may lead to potential errors inpredicting the cell stiffness [170].3.3.2 Effect of the cell geometryIn the ring stage, the malaria-infected RBC retains the shape and size ofthe healthy RBC. Advancing into the trophozoite and the schizont stages,however, it undergoes dramatic morphological changes (cf. Table 3.1). Over-all, the transformation can be viewed as a steady loss of the excess surfacearea, as the cell becomes increasingly swollen [196]. In this subsection, weinvestigate how such geometric changes alter the iRBCs’ deformability. Ac-cording to Nash et al. [153], the shear modulus of the membrane remainslargely constant through the ring, trophozoite, and schizont stages, around15 to 25 µN/m. Here, we fix the shear modulus to be Gs = 15 µN/m, andexamine the passage and blockage of cells having the four shapes indicatedin Table 3.1 for the ring, early trophozoite, late trophozoite and schizontstages. Although we will refer to these cell morphologies by the stages ofinfection, the cells have no merozoite inside.Consider the narrow channel first, with a width of w = 3.2 µm in themid-section. Our simulations show that only the ring-stage iRBC can pass,and all the other cell shapes block the channel. Figure 3.8 depicts the steady-state shape of the iRBC after blockage. Once the iRBC is arrested at thecontraction, the negative pressure downstream draws out a tongue from thecell membrane, whose length lt is plotted as a function of the excess surfacearea ratio se. With the loss of excess surface area, lt decreases monotoni-cally, indicating the gradual loss of the cell’s deformability. Increasing thepressure drop ∆P or capillary number Ca will deform the cell more andproduce a longer protrusion. However, within the range tested in our nu-merical experimentation, we were unable to deform any of the three shapessufficiently so that the cell passes the narrow channel.In this connection, it is interesting to note the concept of minimum cylin-drical diameter (MCD) that has been used to interpret recent experiments353.3. Resultsschizontlate trophozoiteearly trophozoite0. 3.8: Steady-state shape of the iRBC in the early trophozoite, latetrophozoite and schizont stages after it blocks the narrow channel (w =3.2 µm). The length of the tongue lt, defined as the distance between thecell front and the channel entry, and scaled by the channel length L, isplotted as a function of the excess surface area ratio se. Color contours ofthe surface stretching γ are also shown in the insets.on cell entry into microfluidic channels [4, 77]. Because of its fixed surfacearea and volume, a red blood cell can at most be stretched into a cylinderof a certain length and diameter, and thus would fail to enter a thinnercapillary regardless of the magnitude of forcing. In our view, the MCD thusdefined is not a very useful indicator of cell blockage of capillaries. It is apurely geometric lower-bound that does not account for the dynamic aspectsof the passage process. For instance, in microfluidic channels as well as invivo, only finite pressure gradient is available for deforming the cell. Thus,blockage may appear in channels much wider than MCD, far before the cellis stretched into a cylindrical shape (cf. insets in Figs. 3.6 and 3.8). Nu-merical experiments show this to be the case, especially for relatively highvalues of the membrane modulus. Besides, MCD is defined for a cylindricaltube, and does not apply to passage through the narrow slit in the spleennor through microfluidic channels having rectangular cross-sections.In the medium and wide channels, most of the four iRBC shapes pass363.3. ResultsFigure 3.9: The transit time tt through the medium and wide channels asfunctions of the cell excess surface area ratio se. The transit time is scaledby the T .without occlusion. The sole exception is the schizont-stage iRBC blockingthe medium channel. Figure 3.9 plots the cell transit time tt as a functionof se. From the ring stage (se = 1.34) to the late trophozoite stage (se =1.15), the cell surface area remains more or less constant while the cellvolume increases (Table 1). Consequently, tt increases monotonically withdecreasing se in both channels as the gradual swelling of the cell makesdeformation more difficult. The schizont shape, more spherical with thelowest se = 1.04, corresponds to much reduced cell volume and surface area.This decreases deformability and causes blockage of the medium channel.Through the wide channel, on the other hand, not only does the schizont-iRBC manage to pass, but it takes the least transit time among all cases.This anomaly is mainly due to the smaller cell volume of the schizont shape;the size of the channel is such that the schizont iRBC barely needs to deformto enter the contraction. To sum up, the excess surface area ratio se isa key parameter to cell deformability, especially when large deformationis required for passing narrow pores. In relatively wide channels, the sebecomes a secondary factor, while the cell size plays a key role.373.3. Results3.3.3 Effect of the parasites clusterThrough the three stages of malaria infection, the parasite changes shapeand increases in size, eventually taking up more than half of the volume ofthe host red blood cell in the schizont stage [47, 192]. In this subsection, wefocus on the effect of the parasite volume on transit of the iRBC througha microfluidic channel. In the ring stage, the merozoite develops into adisc 2 µm in diameter and 0.5 µm in thickness. Simulations show that thering-stage iRBC can pass through all three channels. This agrees with theobservations of Shelby et al. [195] in microfluidic channels of similar size(Fig. 3.1), and is also consistent with the conclusion of Hosseini and Feng[87] that the ring-stage merozoite is too small to hinder cell deformation instretching simulations.In trophozoite and schizont stages, we have systematically varied theparasite size for each of the iRBC shapes of Figs. 3.3 and 3.4 and Table 3.1,with the membrane shear modulus fixed at Gs = 15 µN/m. A clear trendemerges for all three channels: increasing parasite size eventually causesblockage, and the minimum parasite size increases for blocking wider chan-nels. In the following we will only present results of the schizont iRBC in thewide channel (w = 4.8 µm), with the iRBC having a spheroidal shape andcontaining a rigid spheroidal parasite of the same aspect ratio 1.6 (Fig. 3.4).In reality, the parasite multiplies to dozens of new merozoites that are closelypacked inside a digestive vacuole at the center of the iRBC [71, 210]. Ourspheroidal parasite represents the entire vacuole. The iRBC radius is fixedat 3.2 µm, and we vary the radius of the parasite from 1.6 to 2 µm.Figure 3.10 shows the trajectories of the iRBC for two sizes of the par-asite. A baseline case with no parasite is also shown for comparison. Notsurprisingly, having the parasite slows down the passage. With the largeparasite, the iRBC blocks the contraction. Note that in this case the par-asite is still considerably smaller than the channel size. The occlusion isnot due directly to the solid parasite blocking the channel, but indirectlyto the parasite constraining the deformation of the cell. The membranepresses against the solid parasite and deforms less as a result, a scenariopreviously analyzed in stretching of iRBC by optical tweezers [87]. Recallthat membrane-parasite overlap is prohibited in our model by a repulsivepotential. Not surprisingly, the critical parasite size for occlusion decreasesfor narrower channels and higher membrane stiffness.It is interesting to point out that a smaller parasite tends to move forwardinside the cell during the transit, and ends up in the front part of the cell(see top inset in Fig. 3.10). A larger parasite, on the other hand, feels the383.3. ResultsexitentryFigure 3.10: Trajectories of the schizont-stage iRBC in the wide channel(w = 4.8 µm) with a parasite of different sizes. R = 0 indicates a baselinecase without parasite.resistance from the constricting channel walls. It tends to move more slowlyand shift toward the rear portion of the cell. In particularly, when blockagehappens, the parasite is pushed against the rear end of the membrane (seebottom inset in Fig. 3.10). Conceivably this is a precursor to “pitting”,when the parasite is expelled from the cell through a temporary openingin the membrane. Such a process has been observed in vivo [185] and inmicrofluidic channels [195].3.3.4 CombinationsIn this subsection, we integrate all three factors examined above, stiffen-ing of the cell membrane, changes in cell morphology and parasite growth,into a series of simulations that mimic the evolution of iRBC through the var-ious stages of infection. Following microfluidic experiments [4, 67, 77, 195],we correlate the appearance of blockage in channels of different width withthe overall loss of deformability. For the healthy red cell, we take Gs = 5µN/m, while for all stages of infection, we use a constant Gs = 15 µN/m393.3. ResultsTable 3.2: The passage and blockage of healthy and infected red blood cellsthrough microfluidic channels of different size.w = 3.2 µm 4 µm 4.8 µmHealthy passage passage passageRing passage passage passageEarly trophozoite blockage blockage passageLate trophozoite blockage blockage passageSchizont blockage blockage blockagein accordance to prior measurements and computations [63, 87, 153]. Theshape and size of the iRBC and the parasite are those of Table 3.1 andFigs. 3.2–3.4. Note that we differentiate between the early and late tropho-zoite stages.Inspired by images from the microfluidic experiments of Shelby et al.[195] (Fig. 3.1), we summarize our results in a tabular form (Table 3.2). Ev-idently, with the progression of malaria infection, the iRBC gradually losesits deformability and causes blockage of larger and larger channels. Thehealthy RBC is able to pass through all three channels used in the currentstudy. So is the iRBC in the ring stage. Infected red cells in the earlyand late trophozoite stages pass through the wide channel but block themedium and narrow ones. In the schizont stage, iRBCs block even the widechannel. These results are in qualitative agreement with the experimentalresults of Shelby et al. [195] In fact, one may even claim a degree of quanti-tative agreement. For example, they reported that trophozoite-stage iRBCsblock a channel 4 µm in width but pass freely through a 6 µm channel,and schizont-stage iRBCs block a 6 µm channel but traverse a 8 µm chan-nel. However, we should note that our channel geometry and pressure drophave not been precisely matched with those in the experiment. Incidentally,Shelby et al. [195] also observed an interesting phenomenon of “jamming”,with two or more cells clogging the contraction leading to the narrow partof the flow conduit (Fig. 3.1). In this case, the adhesion among iRBCs playsan important role in causing the blockage. Our simulation deals only withsingle red cell, the jamming phenomenon involves multiple cells interactions,which might be our future effort.403.4. Conclusion3.4 ConclusionIn this chapter, we report a smoothed-particle-hydrodynamics simula-tion of the transit of healthy and malaria-infected red blood cells throughmicrofluidic channels. The cell membrane is represented by a particle-springnetwork that possesses in-plane elasticity as well as bending elasticity. Theparasite is modelled by a group of particles constrained to move as a rigidbody. The cytosol and surrounding liquid are discretized into standard SPHparticles, and the entire fluid-solid system is simulated by an explicit SPHalgorithm. With progression of the malaria infection, the red cell graduallyloses its deformability, measured by its ability to pass through a microchan-nel of a certain dimension. The main objective of the work is to probe howchanges in the cell membrane rigidity, the excess surface area of the cell andthe size of the parasite, individually and taken together, affect the transit ofthe infected red cell through a microfluidic channel. Within the parameterranges explored, the results can be summarized as follows.(a) The rigidity of the cell membrane hinders passage of the cell throughmicrofluidic channels. The transit time increases with increasing mem-brane modulus, especially for the narrower channels. Blockage occursfor sufficiently stiff membranes. These results are in qualitative agree-ment with experimental measurements [67].(b) The progressive loss of excess surface area impairs the capability of theinfected red cell in traversing narrow channels, and blockage occurs forsufficiently low excess surface area ratio. This agrees qualitatively withprevious microfluidic experiments [76, 77].(c) The presence of the parasite reduces the deformability of infected redcells, and the effect increases as the parasite grows. The parasite maycause occlusion of microfluidic channels considerably larger than theparasite size. This result confirms previous speculations [153], and isconsistent with prior computational studies of cell stretching [87].(d) In reality all three factors are at play, and the numerical simulationssuccessfully reproduce the experimental picture of infected red cellsoccluding larger and larger channels with the progression of infection.However, it is not straightforward to quantify the separate contributionof each factor and make a meaningful comparison among them.(e) The numerical simulations suggest that red cell rigidification at dif-ferent stages of infection can be quantified by measuring the critical413.4. Conclusionpressure drop required for passage through a constricting channel ofa certain geometry. This provides a basis for developing microfluidiccell assays.We should point out the assumptions and simplifications of our physicaland numerical model that may have limited the fidelity and accuracy of thenumerical results. The actual RBC membrane has a complex structure, andonly the grossest features are retained in our particle-spring network. Assuch, areal conservation is enforced via an ad hoc energy penalty, not fromthe intrinsic constitutive behavior of the membrane. A potential remedy tothis is to employ nonlinear springs that generate stronger strain-hardeningfor the membrane as a whole [160]. In addition, we have omitted mem-brane viscosity altogether. Prior work has indicated an essential role forthis viscosity in determining the times scale of the membrane relaxation[53]. Moreover, we have assumed equal viscosity for the cytosol and sus-pending fluid for simplicity. In reality, the cytosol viscosity is considerablyhigher than that of the plasma [20], and may evolve dynamically as the par-asite modifies the chemical composition of the cytosol. Finally, this studyhas focused on the deformability of infected red cells and disregarded cy-toadherence, another important factor in the pathogenesis of malaria [38].Pathologically enhanced adherence among infected red cells and with thevascular endothelium plays a definite role in obstruction of microcapillaries,as experiments in microfluidic channels have already suggested [195]. In fu-ture work, these simplifications should be examined and possibly removedin more sophisticated and faithful simulations.42Chapter 4Modelling themechanosensitivity ofneutrophils passing throughmicrofluidic channels4.1 IntroductionNeutrophils are a primary member of the immune system. In its in-active state in vascular flow, the neutrophil has a spherical shape with adiameter of about 10 µm [58, 80, 211]. When the neutrophil gets close toan inflammation site, it will sense chemical signals and become activated.Upon activation, the neutrophil deforms and extends pseudopods. Underthe guidance of the chemical signals, it passes through endothelial slits toreach the inflammation site [104]. In this process, the activated cell under-goes severe mechanical deformation and dramatic morphological changes.Softening of the neutrophil facilitates its passing through the endothelialslits. Besides, the neutrophil activation also supplies the protrusion forcethat drives the extravasation [173]. Thus the biophysical properties of theneutrophil is crucial to its function.Aside from chemical signals, neutrophils are also sensitive to mechanicalforcing and deformation [24, 126, 140, 237, 238, 244]. For example, Cough-lin and Schmid-Scho¨nbein [24], Moazzam et al. [140] and Makino et al. [126]investigated the response of the neutrophils to shear stresses. In their experi-ments, the neutrophils are initially activated due to interactions between themembrane and substrate integrins. That is to say, the neutrophil presentspseudopods and is actively exploring its surroundings. If a mild shear stressis applied by imposing a gentle fluid flow over the neutrophil, the protrusionsretract, and the neutrophil stops crawling and recovers its spherical shape.These studies suggest that a mild shear stress inhibits neutrophil activation.Applying a much stronger mechanical force, Zhelev and Hochmuth [244] used434.1. Introductiona micropipette to aspirate the neutrophil membrane, which leads to severedeformation of the neutrophil and detachment of the cortex network fromthe plasma membrane in the portion of the cell membrane that has beensucked inside. Interestingly, this triggers actin polymerization and stiffeningof this portion, and the neutrophils were activated with formation of protru-sions. This points to stress-induced activation of the neutrophil. Later, Yapand Kamm [237, 238] studied the transit of neutrophils through microflu-idic channels and filters. They found that the mechanical deformation maycause activation of the neutrophils with the formation of pseudopods.These seemingly contradictory responses of neutrophil to mechanical de-formation [24, 126, 140, 237, 238, 244] can be rationalized by the FLNa-FilGAP pathway [46, 197]. Filamin A (FLNa) is a binding protein thatconnects the F-actin filaments, and has a binding site for recruiting theFilGAP proteins from the cytosol. FilGAP is notable as an inhibitor forRac [46, 152, 158, 197], which regulates protrusion of the pseudopodia andlamellipodia during neutrophil migration. Once deformed by external forc-ing, FLNa is known to release FilGAP to the membrane where it inhibitsRac [46, 197]. When neutrophil is under mild forcing, such a FLNa-FilGAPpathway leads to inhibition of Rac. Hence, neutrophil activation is abol-ished [24, 126, 140]. In the micropipette aspiration experiments of Zhelevand Hochmuth [244], on the other hand, the severe cell deformation causesthe actomyosin network to detach from the plasma membrane. This voidsthe FLNa-FilGAP pathway. Thus, Rac remains highly active on the mem-brane, and induces the cortex to grow back and new protrusion to form.This explains the observation that neutrophils activation is inhibited undermild shear stress, but activated under severe deformation.However, the experiments of Yap and Kamm [237, 238] raise furtherquestions on how the neutrophil responds to mechanical forces. In forcingneutrophils through microfluidic channels with a cross-sectional dimensionon the order of a few microns (Fig. 4.1), Yap and Kamm [238] found that theneutrophil fluidizes upon entering the narrow channel. That is, its cortexdisassemble and its elastic modulus drops sharply. This allows the neu-trophil to deform and “flow” into the narrow channel as if a liquid drop.They also found that temperature has little effect on fluidization. This im-plies that fluidization is a purely mechanical process that does not involvechemical rates. Once inside the microfluidic channel, the neutrophil rebuildsits cortex in time, and then develops pseudopods to explore its surround-ings. Thus, the neutrophil is activated. Temperature strongly influences thetime required for the neutrophils to activate, and this is taken to mean thatactivation is regulated by chemical signals such as the Rho GTPases.444.1. IntroductionFigure 4.1: Microfluidic channel. Cited from Yap and Kamm[238].In a later experiment, Yap and Kamm [237] probed the dynamics ofneutrophil activation by varying the time in which the cell was subject tomechanical deformation, and by imaging the density and distribution of thecortex during activation. This was done in a filtration device that allowedneutrophils to pass through the filter at a slower or faster rate. The filteredcells were collected downstream and the concentration of F-actin is recordedthrough microscopy. As F-actin is the major component of the cortex, it isused to mark the evolution of the latter. Their results indicate that undera fast transition rate, the F-actin content drops suddenly and then growsback uniformly in the cell without formation of pseudopods. Under the slowtransition rate, on the other hand, the cell forms pseudopods where F-actindisplays a high concentration locally. Yap and Kamm [237] noted the effectof the entry time te on neutrophil activation. They postulated that when te islarge enough, the cell has time to be activated while deforming and enteringthe filter. In their experiments, te is about 15 s for the slow filtration rateand 2 s for the fast filtration rate. These time scales are comparable to thereciprocal of the kinetic rate of the Rho GTPases [159, 166, 179], which istypically several seconds.The brief literature review above demonstrates the complexity and sub-tlety in how neutrophils sense and react to mechanical forces and defor-mation. Severe deformation may destroy the cytoskeleton and bring aboutfluidization. Loss of the cytoskeleton removes the FLNa-FilGAP pathwaythat may otherwise inhibit the reconstitution of the cortex. The recovery454.2. Neutrophils modelthereafter can take two possible routes: uniform recovery or activation, thelatter distinguished by formation of pseudopods with localized high concen-tration of F-actin. The activation scenario appears consistent with earlierfindings of Zhelev and Hochmuth [244].This work was motivated by the experiments of Yap and Kamm [237,238]. More specifically, we wish to investigate these questions: (a) How dothe Rac and FLNa-FilGAP pathways govern the activation process? (b) Howis activation affected by the competing time scales of cell entry and chemicalreaction? (c) How does the cell cortex recover with or without activation?Our objective is to synthesize hypotheses and ideas from prior experimentsand modelling into a quantitative framework for describing and predictingthe fluidization and activation of neutrophils. This will be realized in abiomechanical model that couples chemical signalling with cell movementand membrane deformation in a fluid environment.4.2 Neutrophils modelTo investigate the questions raised above, we have carried out two-dimensional (2D) numerical simulation of the neutrophil transit througha narrow channel. It is important that we capture not only the overall me-chanics of cell deformation and fluid flow, but also the intracellular remodel-ing that characterizes fluidization and activation. Therefore, our neutrophilmodel accounts for the elasticity and deformation of the cell membrane, dy-namics of the cytoskeleton, and the chemical kinetics involved in signaling.In the following subsections, we will discuss each of those components atlength. As the nucleus and other cytoplasmic organelles play no role in thephenomenon of interest, they will not be explicitly modeled.4.2.1 Mechanics of membrane deformationThe neutrophil membrane consists of two components: the plasma mem-brane and the cytoskeletal network, also known as the cell cortex. Theplasma membrane can experience only very limited strain (5%) before rup-ture [141]. However, there are wrinkles on the neutrophil membrane in itsresting, spherical shape. The excess surface area allows large morphologicalchanges after activation of the neutrophils. The cortex is a thin layer ofactomyosin network underlying the plasma membrane, and consists princi-pally of F-actin filaments and myosin motors. In the rest state, the cortexis pre-stressed due to relative sliding of the myosin motors among the poly-mer filaments. In this study, the cell cortex never detaches from the plasma464.2. Neutrophils modela) b)Figure 4.2: (a) Schematic of the particle-spring model for the cell membrane.The portion inside the red box is shown in the magnified view on the right.(b) Xi denotes the location of the ith particle, li is the length of the ithedge, and Θi is the angle between the ith and i+ 1th edges. At the start ofsimulations, all edges have the same length l0, and all neighbouring edgeshave the same angle Θ0 = 2pi/N .membrane. Thus we need not represent each separately. By the same token,the term “membrane” will refer to both components hereafter, unless eithercomponent is mentioned explicitly.In our 2D model, schematically depicted in Fig. 4.2, we represent themembrane by N linearly elastic springs connected at N nodes into a loop.In the initial, undeformed shape, the nodes fall on a circle of diameter d.The plasma membrane contributes in-plane tension as well as out-of-planebending forces. The in-plane tension on the ith edge is thus given as:σim = Em(li − l0l0), (4.1)where Em is Young’s modulus, li is the instantaneous length of the ith spring,l0 is the resting length of all springs. Note that in our 2D representation, thein-plane stress consists only in a tension. While the cell membrane is oftenassumed to have nonlinear elasticity [39, 128], we use the linear strain-stressrelationship for simplicity. For the out-of-plane bending force, we assign thefollowing bending energy[216, 224, 235]:Eb =kb2∑j=1,Ntan2(Θj2), (4.2)474.2. Neutrophils modelwhere kb is the bending modulus. Thus, we assume zero spontaneous cur-vature in the membrane. The derivative of the energy with respect to nodalposition gives rise to the bending force (see Eq. 4.6 below).Cortical tension is produced by myosin dynamics. This is typicallytreated as a tension force on the elastic membrane in numerical studies[229]. Thus, we assume that the cortex contributes an additional tensionbut no bending. For the ith edge, this cortical tension is given byσiτ = Eicεi, (4.3)where Eic is the elastic modulus for the ith cortical segment, and εi is thepre-strain of the cortex due to myosin motors [45, 162]. Note that Eic isnot constant but varies along the cell membrane and in time during cellfluidization and activation. It will be used to indicate the integrity of theactomyosin cortex. The myosin phosphorylation is activated by the myosinlight chain kinase (MLCK), which is regulated through the small GTPasesRhoA pathway [103], and thus the pre-strain εi is related to the local Raclevel in our polarization model. More details will be given in the nextsubsection.We also apply a volume constrain to ensure that the cell volume, whichis really the cell area in our 2D model, is conserved to within 1% during thesimulation. This is implemented through an energetic penalty:Ev = kv(V − V0V0)2. (4.4)This constraint is motivated more by numerical than physical considera-tions. In principle, incompressibility for the fluid inside and outside the cellshould imply cell-volume conservation. Because of numerical errors, how-ever, volume variation can become a concern when the cell undergoes severedeformation. Hence, this “superfluous” constrain has been used in red bloodcell models in the past[36, 55, 235], and will be applied to our neutrophilmodel as well.Lastly, we include a protrusion force on the cell membrane that is dueto the growing cytoskeleton underneath. This is intended to capture theformation of pseudopods during activation. Following previous models [218],we express the protrusion force as a function of the active form of the signalprotein Rac on the membrane:F ipro = kproAia0l0ni −kproN∑j=1,NAja0l0nj , (4.5)484.2. Neutrophils modelwhere kpro is a protrusive coefficient, Ai is defined by 12(aili + ai+1li+1),where ai is the concentration of Rac on the ith edge and li is the length ofthe ith edge. a0 is the initial concentration of Rac on the membrane. ni isthe unit outward normal on the ith node, which bisects the angle betweenthe adjacent edges. The first term is the protrusion force intended for theforefront of the cell, where Rac is expected to build up (more details insubsection 4.2.2 below). The second term represents the reaction force. Asthe neutrophil will not be anchored onto a substrate via focal adhesion, theprotrusion force must be balanced by backward forces on the membrane inthe rear. Thus, the sum of F ipro over the cell periphery vanishes.Collecting all the forces discussed above, we write the total force actingon the ith membrane node asF i = F im + Fiτ + Fipro +∂(Eb + Ev)∂Xi. (4.6)where F im = (σi+1m τ i+1−σimτ i) and Fiτ = (σi+1τ τ i+1−σiττ i) are the tensionson node i due to the plasma membrane and cortex, respectively. τ i denotesthe tangent vector pointing from the (i− 1)th point to the ith point.4.2.2 Kinetics of chemical signalingTypically, the polarity of the neutrophil is characterized by anisotropicdistribution of the Rho family of GTPases, such as Rac, RhoA and Cdc42.Prior experimental data and modelling show an antagonism between Racand RhoA that is modulated by Cdc42, and that Rac serves as an indicatorof cell polarization [119]. Thus, we will use Rac as the sole marker forneutrophil polarity in our model, and will not account for RhoA and Cdc42separately.In the model, we assume that the chemical signal Rac has two states: anactivated state that binds to the membrane, and an inactivated state thatfloats in the cytoplasm. The active Rac is defined as a concentration ai oneach edge of the membrane (Fig. 4.2b), while the inactive cytosolic Rac isgiven by I. The two are related by the conservation of the total amount ofRac:I +N∑i=1aili = 1. (4.7)Thus, both ai and I have been scaled by the total amount of Rac, such thatI is dimensionless, and ai has the dimension of the reciprocal of length.494.2. Neutrophils modelThe active Rac evolves according to a kind of reaction-diffusion equation:d(aili)dt= k+lidI − k−(aili) + (qi − qi−1), (4.8)which is essentially a “mass balance” equation for the total amount (aili) ofactive Rac on the ith edge. The first term on the right hand side representsattachment and activation of inactive Rac I onto the edge; d is the initialcell diameter and (li/d)I gives the inactive Rac available to the edge. Therate of Rac activation is written ask+ = kon +γa2iK2 + a2i, (4.9)where the Hill function is introduced to stabilize a polarized state by wave-pinning [127, 146, 147]. The second term of Eq. (4.8) represents down-regulation of the membrane-bound active Rac, withk− = koff + kτσiτ (4.10)incorporating the inhibition of Rac by cortical tension στ through the FLNa–FilGAP pathway [46, 197]. The last term on the right hand side of Eq. (4.8)represents 1D diffusion of active Rac along the membrane, withqi = Dai+1 − ai12 (li + li+1)(4.11)being the diffusive flux from edge li+1 to edge li at node i, D being thediffusivity. As the signaling proteins diffuse about 100 times faster in thecytosol than on the membrane [166], we do not account for diffusion ofI in the cytosol but assume that the inactive I is uniformly available toattachment onto the membrane and activation. With the above, Eq. (4.8)can be written as an evolution equation for ai(t):a˙i =(kon +γa2iK2 + a2i)Id− (koff + kτσiτ )ai −ai l˙ili+D∇2i ai, (4.12)where the dot indicates time derivative, and ∇2i is the discretized form ofthe Laplacian:∇2i ai = 2ai+1(li + li−1)− ai+1(li−1 + 2li + li+1) + ai−1(li + li+1)(li + li+1)li(li + li−1). (4.13)Once the neutrophil is polarized, it is known that Rac is enriched inthe front and promotes F-actin growth, while RhoA is concentrated in the504.2. Neutrophils modelrear and promotes myosin contraction there [234]. Since our model does notexplicitly account for RhoA or myosin, we represent the polarity in myosincontraction indirectly via the pre-strain εi (Eq. 4.3):εi = ε(ai) = ε0 exp[1−aia0i]. (4.14)Thus, the pre-strain ε, and the cortical tension στ by extension, are inhibitedby the high Rac in the front, and promoted by the low Rac (and implicitlyhigh RhoA) at the rear. The pre-strain of the resting state ε0 is taken to bea constant.4.2.3 Remodeling of the cytoskeletonAs mentioned above, our model does not represent the cytoskeleton bya physical element separate from the plasma membrane, but uses the cor-tical modulus Ec to indicate the density and integrity of the cortex. Thedynamics of Ec is regulated by mechanical perturbations as well as chemicalsignals (i.e., Rho family of GTPases). It can be reduced to zero by severedeformation, and then grow back under the promotion of Rac.Alvarado et al. [3] reported that mechanical forces acting on the corticalnetwork causes rupture of the bonds between the actin filaments. Costaet al. [23], Murrell and Gardel [151] and Sato et al. [182, 183] reported thatthe actomyosin network disassembles much more easily under bending thanstretching. In our model, therefore, we assume that the cytoskeleton dis-assembles catastrophically when the instantaneous bending rate exceeds athreshold:Ifθ˙iΘ0> β, then Eic = 0, (4.15)where β is a critical bending rate constant for disassembling of the network,and Θ0 is the initial angle between adjacent segments. θ˙i is the bending rateon the ith edge, which is defined from the bending rate of the end nodes asθ˙i = 12(|Θ˙i−1|+ |Θ˙i|). This is how our model realizes fluidization of theneutrophil.A primary role of active Rac on the membrane is to induce actin poly-merization and strengthen the local actomyosin network [124]. Furthermore,Yap and Kamm [237] documented localized high concentration of F-actinat the protrusion front, which is consistent with the idea that active Racpromotes F-actin polymerization and growth. Based on these results, wemodel the recovery of the cortex after fluidization toward the homeostatic514.3. Numerical methodslevel byE˙ic = kpolyaia0(E0c − Eic), (4.16)where kpoly is a polymerization rate for the cytoskeletal network, ai and a0are instantaneous and initial concentration of the signal protein Rac.4.3 Numerical methodsThe numerical computation consists of two major elements: the fluidflow inside and outside the neutrophil, and the deformation of the cell mem-brane. In a mechanical context, this is essentially a fluid-structure interac-tion problem. These two elements are supplemented by the evolution of Racand cortical modulus Ec on the cell membrane, which is a one-dimensionalproblem and thus numerically simpler.4.3.1 Fluid dynamicsThe fluid dynamics is governed by the incompressible Navier-Stokesequations:∂v∂t+ v · ∇v = −1ρ∇p+µρ∇2v + fw + fm, (4.17)∇·v = 0, (4.18)where v is the velocity, ρ is the density, and µ is the viscosity of the fluid.fw and fm are forces acting on the Eulerian nodes due to the immersedsolid wall and membrane, respectively. These will be discussed next.4.3.2 Immersed boundariesThis study involves both elastic and rigid boundaries in contact withfluid. We treat the fluid-solid interaction by the immersed boundary method[52, 163]. The Lagrangian points (Xi) do not necessarily coincide withthe Cartesian nodes (x), of course. Thus, a smoothed delta function isapplied to distribute the elastic membrane forces to the surrounding fluidmedium. For the stationary wall boundaries, the “direct forcing” algorithmis applied to enforce the fluid motion to satisfy the non-slipping condition.For this scheme, the target velocity on a Eulerian node that close to thewall boundaries is calculated, then the “direct-forcing” on this node will bedirectly computed. Details have been given in Chapter 2, Methodology.524.4. Simulation setup and model parametersFigure 4.3: Geometry for the channel with a constricting section.We have validated the flow solver by computing the deformation of acapsule in a shear flow, and comparing the results with those in the literature[203]. We refined the Eulerian grid size h as well as the membrane resolutionl0, with the following 3 test cases: (i) h = 1/32 and l0 = 1.4h, (ii) h = 1/64and l0 = 2.8h and (iii) h = 1/64 and l0 = 1.4h. The results of all 3 casesagree with those of Sui et al. [203] within 2%. The following microfluidicsimulations use the mesh scheme (iii).4.4 Simulation setup and model parametersThe computational domain for our simulations is depicted in Fig. 4.3.It is a rectangular channel of dimensions Lx × Ly, with a constriction. Theneutrophil initially has a circular shape with diameter d. The constrict-ing channel has width w and length Lc. Periodic boundary conditions areapplied in the flow direction, between the entrance and the exit, whereasno-slip conditions prevail on the solid walls. A pressure drop ∆P is imposedacross the length of the domain Lx, which sustains the channel flow. Wetake the characteristic length to be L0 = Ly, and define a characteristicvelocity by U0 = w28η∆PLx. The characteristic time is then T0 = L0/U0.The model parameters are summarized in Table 4.1, with literaturesources given where available. They fall into 3 groups: geometric, physi-cal and kinetic. In the following, we briefly discuss the choice of parametervalues in the simulations.For the geometric parameters, we set Lx = 80 µm and Ly = 20 µm. Inthe filtration system of Yap and Kamm [237], the pore diameter ranges from3 to 5 µm, and the pore length varies from 9 to 10 µm. In our simulation,we choose w = 5 µm and Lc = 10 µm, respectively. The length of the widechannel upstream of the constriction is 2Lc, and that downstream is 5Lc.534.4. Simulation setup and model parametersTable 4.1: Parameters for the model. The value for the kpoly is calibratedthrough comparing the simulation with the experiment (see discussion inSubsection 4.5.3).Symbol Meaning Value and referencesLx length of the channel 80 (µm)Ly width of the channel 20 (µm)w width of the constricting section 5 (µm) [237, 238]Lc length of the constricting section 10 (µm ) [237]d diameter of the cell 10 (µm ) [58, 80, 211]η cytosol & fluid viscosity 0.5 (Pa·s)ρ cytosol & fluid density 1000 (kg/m3)∆P driven pressure drop 25 ∼ 150 (Pa) [237, 238]σ0 cortex tension 50 (µN/m) [73, 80, 244]ε0 pre-strain 0.1 [23, 121]Em membrane modulus 500 (µN/m)kb bending coefficient 1×10−15 (J)β cortex network disassembly criterion 10 (s−1)kpro protrusion coefficient 8 (pN) [169, 218]D diffusivity 0.5 (µm2/s) [26, 166]kpoly polymerization rate 0.02 (s−1) [237]γ wave-pinning parameter 0.0125∼0.375 (s−1)K wave-pinning parameter 0.0187 (µm−1)The initial diameter of the neutrophils is set to be d = 10 µm.Among the physical parameters, the viscosity ratio between the cell cy-toplasm and the surrounding fluid is set to be 1, and the viscosity is chosento be η = 0.5 Pa·s. In reality, the cytoplasmic viscosity of neutrophils ismore than 104 times greater than that of the suspending fluid [214]. Sim-ulating such a large viscosity ratio is numerically challenging, if possible atall [171]. Thus we have adopted the artificial viscosity value that is betweenthose of the cytoplasm and the suspending fluid, such that for the range ofpressure drop across the microfluidic channel, the neutrophil will enter andpass through the microfluidic channel on time scales comparable to thosein the experiments of Yap and Kamm [237, 238]. The pressure drop ap-plied to push the cell through the constricting channel will be varied withinthe range ∆P = 50 ∼ 200 Pa after the experimental conditions of Yap andKamm [237, 238]. The Reynolds number Re = ρL0U0/η is vanishingly smallfor the above parameters (∼ 10−7), but the inertial term is retained in the544.5. Results and discussionmomentum equation.The measured values of the rest-state cortical tension σ0 falls between30 and 90 µN/m [73, 80, 244]. Hence we set σ0 = 50 µN/m for the neu-trophil at resting state. The prestrain of the resting state ε0 is chosen to be0.1 according to the literature [23, 121]. Now the cortical modulus in theresting state can be backed out: E0c = σ0/ε0 = 500 µN/m. Note that inthe cortical tension the plasma membrane contribution has been excluded.In the study of Yap and Kamm [238], the residual elastic modulus after theneutrophil fluidization is about 50% of the original one, which is mainly dueto the plasma membrane. Thus we set the elastic modulus for the plasmmembrane to be the same as that of the cortex Em = 500 µN/m. The majordimensionless parameter that controls the cellular dynamics is the capillarynumber Ca = ηU0/σ0 defined using the rest-state cortical tension. Usingparameters in Table 4.1, we have Ca = 0.0195 ∼ 0.078 in our simulations.The value of the β is chosen such that the small deformation of the cellin the flow does not disrupt the cytoskeleton, but large deformation dueto the microfluidic channel is able to cause the cytoskeleton to disassem-ble. The protrusion coefficient kpro is chosen after a previous model [218].With this value, the maximum value of the protrusion force calculated bykpro(amax/a0) is about 180 pN/µm in our simulation, which matches exper-imental data of Prass et al. [169] (∼ 100 polymerizing actin filaments permicrometer of the leading edge with each filament producing roughly 2 ∼ 7pN of force).The diffusion rate D of active Rac on the membrane ranges from 0.1 to0.5 µm2/s [26, 166], and we use D = 0.5 µm2/s in our simulation. In thesimulation, we will explore a range of γ values while keeping the followingfixed ratios of the coefficients: kon : γ : koff : kτσ0 = 0.05 : 1 : 2 : 0.4. Thisstrategy is informed by the fact that Eq. (4.12), if divided by γ on bothsides, will contain the ratios of these kinetic rate constants on the right-handside, with γ scaling the kinetic time scale on the left-hand side. K is thesaturation parameter [146, 147], and its value is chosen to allow the systembistable solutions and the wave-pinning mechanism for polarization. Basedon these parameters, the equilibrium Rac concentration on the membranebefore cell deformation is a0 = 2.67× 10−3 1/µm.4.5 Results and discussionIn their experiments, Yap and Kamm [237, 238] observed two majorepisodes during the neutrophil’s passage through a narrow channel. One is554.5. Results and discussiona sudden drop of the cell modulus when the neutrophil enters the channel,which is called fluidization. The other is the activation of the neutrophil withformation of pseudopods. Pseudopod protrusion is necessarily preceded bypolarization of the key signal proteins, especially Rac and RhoA. Accord-ingly, we divide our presentation below into three subsections on fluidization,polarization, and activation of the neutrophil. Note that polarization indi-cates a non-uniform distribution of the signal content on the cell membrane,and activation is characterized by a protrusion force pushing on the cellmembrane and causing pseudopods. The latter is a result of the former, asconcentrated Rac on the membrane promotes polymerization of the F-actinfilaments.In our simulation, we fix the ratio among the kinetic rate constants, anduse one constant γ to represent the chemical kinetic rates. We have testedtwo γ values, γ = 2 s−1 for slow kinetics and γ = 4 s−1 for fast kinetics.In the experiments of Yap and Kamm [237], after filtration the neutrophilswere collected in a reservoir without flows. Thus, in the simulation we setthe pressure drop ∆P to zero once the cell’s centroid reaches 1.8d outsidethe exit of the constriction. Then it may relax in an essentially quiescentmedium, as in the experiment.4.5.1 Fluidization of the neutrophilsIn the experiments of Yap and Kamm [237, 238], mechanical deformationcauses a sudden fluidization of the neutrophils due to disassembly of thecortex. In this section, we investigate how mechanical deformation causesfluidization of the cell.Figure 4.4 illustrates the transit of the cell through the channel. Thecolors on the membrane indicate the the cortex modulus Ec, scaled by theresting value E0c . As the cell front deforms and enters the channel, thecortex starts to disassemble at the front first. After the cell has fully enteredthe channel, the cortex has disassembled completely. A similar scenario offluidization occurs for γ = 2 s−1, and so we do not show the results here.This is consistent with the experimental observation that the kinetic rateshave little effect on the fluidization of the cell or the mechanics of its transit[238]. In our model, this is expected of course, as the criterion for corticaldisassembly is based on the instantaneous rate of bending (Eq. 4.15).564.5. Results and discussiona) b)c) d)e) f)Figure 4.4: Passage of the neutrophil through the channel for γ = 4 s−1 at(a) t = 1.54 s, (b) t = 2.30 s, (c) t = 3.07 s, (d) t = 3.84 s, (e) t = 4.61 sand (f) t = 5.38 s. The colors represent the cortex modulus.4.5.2 Polarization of the neutrophilsYap and Kamm [238] found that neutrophils form pseudopods at one endof the cell after the passage through a narrow microfluidic channel, indicat-ing activation. Later, Yap and Kamm [237] observed that the neutrophilswere activated under a slow filtration rate but not under a faster filtrationrate. They hypothesized that the distinction might be because the slow fil-tration rate allows the cell more time to develop polarity in the signallingproteins. It is well known that formation of pseudopods is due to the F-actin574.5. Results and discussionpolymerization, which is regulated by the Rho GTPases Rac [70]. Thus, weneed to investigate the polarization of signalling proteins before neutrophilactivation.It seems necessary to clarify the meaning of polarity in the current con-text. In the previous subsection, we have shown that the cortex fluidizationproceeds from the front toward the end. Where fluidization relieves the cor-tical tension that inhibits Rac on the membrane (cf. Eq. 4.12), the level ofactive Rac a will rise. Hence, the front of the cell has more time for Racgrowth than the rear. Meanwhile, higher a enables the cell cortex (or Ec)to grow back, and the cortical tension will start to inhibit a. Because ofsuch complex interactions, when the cell exits the channel or filter, it nat-urally bears a distribution of a that is higher at one end than the other.Such a spatial gradient is a precursor to polarity. If it is small, it woulddie out quickly by diffusion before activation occurs. If the fore-aft dif-ference is large enough, it can be sustained by a wave-pinning mechanism[83, 127, 130, 146, 147] and persists after the external forcing on the cell hasbeen removed. This corresponds to what has been documented as polariza-tion in the experiments.Figure 4.5 shows the evolution of the Rac polarity for γ = 4 s−1. Beforethe cell enters the channel (Fig. 4.5a), a has a uniform distribution on themembrane: a0 = 2.67 × 10−3 1/µm. Once the cell enters the channel,the deformation initially causes disassembly of the cortex network at thecell front, where a starts to accumulate thanks to removal of the tensioninhibition (Fig. 4.5c). Once a at the cell front surpasses a critical value, thewave-pinning mechanism can amplify a to a high value a+ and maintain thepolarity, as is shown in Fig. 4.5(f). Subsequently, the polarity of a leads to alarge protrusion force at the cell front. Note that in Fig. 4.5(c) and (e), the alevel at the front tip is in fact below the two “shoulder” regions on its sides.This is because the shoulders encountered the corner of the channel first(Fig. 4.5b), where the severe deformation caused the cortex to melt. Thus,a has had the longest time to grow in these regions. Later the shouldersand the front approach a more or less uniform high level (Fig. 4.5f). In thefollowing analysis, we will use the front tip to represent the high-a level inthe front, ignoring the fact that a is maximum at the shoulders.To analyze the polarization process quantitatively, we plot in Fig. 4.6the temporal evolution of the total activated signal A on the membrane,A =∑i=1,Naili, for γ = 2 s−1 and γ = 4 s−1. A starts to increase at t = 1.54s for both γ values. That is when the front of the cell enters the narrow poreand starts to fluidize. For the slow kinetic rate, A reaches the maximum584.5. Results and discussiona) b)c) d)e) f)Figure 4.5: Passage of the neutrophil through the channel for γ = 4 s−1 at(a) t = 1.54 s, (b) t = 2.30 s, (c) t = 3.07 s, (d) t = 3.84 s, (e) t = 4.61s and (f) t = 33.28 s. The colors indicate the level of active Rac a on themembrane.value (A = 0.15) at about t = 5.3 s and then decays slowly to its initialvalue over a time scale of 160 s. The slow, long-term relaxation is not shownin the plot. For the fast kinetic rate, on the other hand, A rises sharplyfirst, peaks at t = 15 s and then stabilizes to a high level: A = 0.29. Thetotal amount of activated signal A on the membrane is higher for the fasterkinetic rate, and it persists without decay for the duration of the simulation,200 s. The cause for this distinction between the slow and fast kinetic rates594.5. Results and discussionFigure 4.6: Evolution of the total activated signal for γ = 2 s−1 (solid line)and γ = 4 s−1 (dash line).will be analyzed in the following.Let us denote the value of a at the cell’s front tip by af and that at theback tip by ab. Figure 4.7(a) shows the evolution of af (solid line) and ab(dash line) for the slow kinetic rate γ = 2 s−1. Before any part of the cortexhas fluidized (t < 1.7 s), a exhibits a uniform distribution, with af = ab.Afterwards, the onset of cortical fluidization at the cell front removes thecortical tension there and allows af to rise. This trend continues until t = 4.3s, when af starts to decline. This reflects the feedback between Rac and thecytoskeleton. The higher af induces cortical reconstitution, which producesa cortical tension that inhibits af through the FLNa-FilGAP pathway. Theactive Rac at the back end of the cell ab follows qualitatively the samebehavior, but with some quantitative differences. First, its rise is delayedby some 3 s relative the the af dynamics. This is because the rear fluidizeslater than the front. Second, when af starts to rise, the cytosolic Rac levelI drops, and this causes a mild dip in ab. Third, as ab starts to rise, thatreduces I and contributes to the rapid decline of af starting from t = 5 s.Finally, ab peaks at the much later time of t = 15 s, and its subsequentdecline is also more gradual. Nevertheless, both ab and af relax to theequilibrium value in about 150 s; the long-term relaxation is not shown in604.5. Results and discussion(a)(b)Figure 4.7: Evolution of the signal at the cell’s front tip (af , solid line) andback tip (ab, dashed line) for (a) γ = 2 s−1 and (b) γ = 4 s−1.Fig. 4.7(a). Thus, this corresponds to a non-polarization case.In contrast, the high kinetic rate γ = 4 s−1 produces strong polarization(Fig. 4.7b). Because of the fast chemical rates, a rises sharply at the shoul-614.5. Results and discussionFigure 4.8: The null-cline for the simplified reaction-diffusion equation(Eq. 4.19). Trajectory (i) corresponds to the polarizing case (γ = 4 s−1)while trajectory (ii) to the non-polarizing case (γ = 2 s−1).ders first (not shown in the plot), and depletes I quickly. Thus, the initialrises of af and ab are both stunted. In particular, by the time the rear endfluidizes (t ≈ 5 s), there remains little inactive Rac in the cytosol to feedthe growth of ab. Thus ab stays low. However, af eventually rises thanks todiffusion from the shoulders. Subsequently a wave-pinning mechanism keepsaf at the high level and ab at the low level. This can be easily appreciatedby analyzing the steady-state solution of the kinetic equation governing a(Eq. 4.12).To illustrate the wave-pinning mechanism, we may ignore for the momentthe diffusion and inhibition of Rac by cortical tension in Eq. (4.12). Thesteady-state solution is thus given by this simplified equation:(konγ+a2K2 + a2)Id−koffγa = 0, (4.19)where a is the local Rac level on the membrane, and I is taken for themoment to be a free parameter, unrelated to a. The steady-state solutiona(I) exhibits the prototypical shape of Fig. 4.8 that allows multiple solutions.Note that the null-cline is the same for the kinetic rates γ = 2 s−1 and 4 s−1,624.5. Results and discussionsince the ratios kon/γ and koff/γ are kept constant in our computations.If I is too small or too large, a single solution prevails, with a low a− orhigh a+, respectively. For the intermediate I values, however, two stablebranches of the solution exist. If a perturbation pushes a past the thresholdmarked by the unstable solution, it will jump from the low branch to thehigher branch (arrow from point 1) or vice versa (point 2). In the actualcomputation, I decreases when a increases thanks to the conservation ofRac (Eq. 4.7). Thus, the trajectory of the solution starts from a large I andsmall a and goes toward the upper left. If a increases sufficiently to crossthe threshold, the solution proceeds to land on the upper branch (trajectoryi). Otherwise it stays close to the lower branch and eventually loops back asthe neutrophil exits the channel (trajectory ii). Both af and ab of the non-polarizing case (Fig. 4.7a) and ab of the polarizing case (Fig. 4.7b) followtrajectory (ii), while af of the polarizing case follows trajectory (i).In the actual simulation, a varies spatially, and its diffusion and inhi-bition by the changing cortical tension modify the simple picture above.The change of the cell circumference during the simulation also plays a role.But the basic mechanism remains the same. In particular, as the corticalfluidization propagates from the front toward the rear, so does the rise ofthe local a. Thus, the a solution exhibits a traveling wave pattern, withthe high a+ solution extending backward in time. This continues until Ibecomes too low to keep raising a to the upper branch. The travelling wavebecomes pinned, and a more or less steady polarized pattern obtains, withhigh a+ at one end and low a− at the other. This is the essence of thewave-pinning mechanism [146].It turns out that the high a+ solution does not necessarily occur at thefront. For certain parameter values, an “inverted polarization pattern” canappear with high a+ at the rear and low a− at the front. This is illustratedin Fig. 4.9 for γ = 3 s−1. As in Fig. 4.7a, af rises first as the front fluidizes,but then declines for the combined effects of cortical tension inhibition andthe rise of ab. Similar to Fig. 4.7a, here ab also surpasses af around t = 7 s.Different from Fig. 4.7a, however, here the chemical rates are high enoughto push ab past the threshold onto the upper branch a+. Hence the invertedpolarization pattern develops.Polarization is not only affected by the kinetic rate γ but also by the entrytime te of the neutrophil through the channel. The latter determines howmuch time is available for Rac to polarize, and thus its role is not unexpected.To vary te, we have simulated four different ∆P values: 50, 100, 150 and200 Pa, which produce the following transit times: te = 4.40, 2.05, 1.31 and0.973 s. These are based roughly on the experimental conditions. In Yap634.5. Results and discussionFigure 4.9: Evolution of the signal content at the cell front (solid line) andback (dash line) for γ = 3 s−1.and Kamm [238], the neutrophil transit time ranges from 0.1 to 10 sec. InYap and Kamm [237], the neutrophil passes through the filters in 2 and 15 sat the fast and slow filtration rates, respectively. This allows us to constructthe “phase diagrams” of Fig. 4.10, with two regimes: no polarization for slowkinetics and/or fast flow, and polarization for fast kinetics and/or slow flow.The fact that the critical γ required for polarization increases with flowrate suggests a criterion based on γte, which may be viewed as the amountof chemical conversion within the cell entry time. Thus, we can recast thephase diagram into Fig. 4.10(b). Although quantitatively the boundarybetween polarization and non-polarization is not horizontal, qualitativelythis supports the argument of Yap and Kamm [237] that polarization isthe outcome of the competition between two time scales. If the chemicaltime γ−1 is within the mechanical dynamic time te, Rac has sufficient timeto develop on the cell membrane, and polarization occurs. Conversely, iftransition rate 1/te is much beyond γ−1, polarization does not have time todevelop before the cell exits the narrow channel.It is tempting to check the above criterion against experimental data.Unfortunately, our γ parameter appears in a highly specific algebraic formfor the Rac evolution equation (Eq. 4.12), and there is no experimental644.5. Results and discussion(a)(b)Figure 4.10: Phase diagram for three regimes of polarization, in the param-eter space of (a) ∆P versus te, and (b) ∆P versus γte.654.5. Results and discussionmeasurement of it. Nevertheless, a rough estimation can be made. Sakoet al. [179] measured the life-time of individual spots of fluorescent Rac onthe cell membrane, and found that more than 80% of the spots disappearwithin 1 s. This suggests a Rac turnover rate of about kRac = 1 s−1. Inthe experiment of Yap and Kamm [237], an entry time of te = 2 s (fastfiltration) does not cause polarization but te = 15 s (slow filtration) does. Acritical value of kRacte between 2 and 15 is consistent with Fig. 4.10(b), insupport of the argument that polarization is the outcome of the competitionbetween the reaction time and the residence time [237].4.5.3 Activation of the neutrophilsYap and Kamm [237] found that the activated neutrophils form pseu-dopods, where F-actin is highly concentrated (Fig. 4.11). It is well knownthat the pseudopod is a result of the protrusion force due to polymerizationof the F-actin filaments, which is modulated by the Rho GTPases includingRac [70]. Experimental images show the pseudopods as bumps on the celledge, whose length scale is considerably smaller than the cell size [237].In the last subsection, we have shown that fast kinetic rates and/or slowtransit of the neutrophil through microchannels may lead to polarization ofthe signal Rac on the membrane, with a high a+ level at the front of the celland a low a− level at the rear. In Eq. (4.5), we have assumed a protrusionforce proportional to the local concentration of a on the membrane. Thus,the main aim of this subsection is to examine the deformation of the cellmembrane under such a force, and the configuration of the resulting pseu-dopod. For this purpose we define a “deformation index” (DI) to reflect thedeviation of cell shape from its original circule shape:DI =√∑i=1,N (2ri − d)2Nd2, (4.20)where ri is the distance of the ith node from the centroid of the cell, d isthe diameter of the cell at the resting state, and N is the total number ofmembrane nodes.Due to the protrusion force acting on the membrane, the polarized cellhas a non-zero DI. Figure 4.12 shows the evolution of the DI at the twokinetic rates examined before, γ = 2 s−1 and γ = 4 s−1. In both cases, DIinitially rises sharply and reaches a maximum of about 0.4 at t = 4.6 s. Thiscorresponds to the point of maximum deformation of the neutrophil in themicrofluidic channel. After that, DI declines as the cell exits the channel664.5. Results and discussionFigure 4.11: Distribution of F-actin inside activated neutrophils afterpassing through filters. The cells are stained with tetramethylrhodamineisothiocyanate-phalloidin, and exhibit highly localized concentrations of F-actin at sites of pseudopod protrusion. (a) shows a neutrophil filteredthrough 3-µm pores under a constant pressure drop, while (b) show a neu-trophil filtered through 3-µm pores under a constant flow rate. Adaptedfrom Yap and Kamm [237] with permission; c©the American PhysiologicalSociety.and starts to retract under membrane elasticity. For the slow kinetic rate,the Rac signal does not polarize and the protrusion force relaxes toward auniform distribution. Thus, DI eventually drops to a value close to zero att = 20 s. This indicates recovery of the cell to its initial circular shape. Forthe fast kinetic rate, the signal polarizes and produces a concentrated pro-trusion force on the front end of the cell (cf. Fig. 4.5c). This is balanced bya much weaker force that is spread over the rest of the membrane. The pro-trusion force sustains a finite cell deformation with DI = 0.07 against elasticrelaxation. The polarized Rac distribution and protrusion will persist, withno appreciable decay within 200 s (Fig. 4.7b). This is our equivalent of aprotruding pseudopod.The distinct long-term behavior of the neutrophil is also reflected by thestate of the cortex, indicated by the cortical modulus Ec. To examine thetemporal evolution of the cortex, we define an average cortical modulusEave =1N∑i=1,NEic. (4.21)674.5. Results and discussionFigure 4.12: Evolution of the cell deformation index (DI) for the slow kineticrate γ = 2 s−1 and the fast rate γ = 4 s−1. The snapshots (a)-(d) show thecell shapes at four interesting times.Figure 4.13 compares the evolution of Eave, together with that of the modu-lus at the front and rear tips, between the activated and inactive cases. Forthe inactive case (γ = 2 s−1), the cortical elasticity grows back more or lessuniformly thanks to the even distribution of the active Rac on the mem-brane. The recovery takes some 140 s to complete. Experimentally, Yapand Kamm [237] reported that for the inactive case the F-actin recovers toits initial level in about 120 s. This apparent agreement of the recovery timeis achieved by tuning the polymerization rate to kpoly = 0.02 s−1. For theactivated case (γ = 4 s−1), the cortex elasticity grows back non-uniformly;it rises much faster at the cell front than at the cell back. This is due to thepolarized Rac distribution in the activated case.Finally, we may compare the activated neutrophil predicted by the modelwith that of experimental observation [237]. After filtration, the activatedneutrophil of Fig. 4.11 shows a non-uniform distribution of F-actin alongthe cell membrane. The maximum F-actin concentration occurs at pseu-dopods. The number and shape of the pseudopods seem quite variable.In Fig. 4.11(a), multiple irregular pseudopods form with dimensions much684.5. Results and discussiona)b)Figure 4.13: Recovery of the cortex elasticity for (a) the inactive case atslow kinetic rate γ = 2 s−1, and (b) the activated case at faster kinetic rateγ = 4 s−1. In each case, we plot the temporal evolution of the corticalmodulus at the front and rear tips of the cell and the average over the entiremembrane. All the moduli are normalized by the resting state modulus E0c .smaller than that of the cell. Fig. 4.11(b), however, features a more reg-ular protrusion on the right side of the cell, together with a more evenlydistributed actin cortex over the top. In Fig. 4.14, we plot our simulation694.5. Results and discussiona) b)Figure 4.14: Predicted geometry and polarity of the activated neutrophilfor γ = 4 s−1 at t = 33.3 s. (a) Distribution of the cortical modulus on themembrane. (b) Distributions of the Rac concentration a and the protrusionforce Fpro indicated by the arrows.results for the activated case. The polarized cortical modulus in Fig. 4.14(a)resembles the polarized actin distribution in the experiment (Fig. 4.11b).However, if we take the front portion of the membrane bearing large a and704.6. ConclusionFpro (Fig. 4.14b) to be the “pseudopod”, then it is much broader and lesspointed than in the experiment. Possibly, the simple 2D representation ofthe cell membrane (Fig. 4.2) and the assumption of a spatially uniform Ihave deprived the model of its ability to capture smaller-scale spatial vari-ations on the membrane.4.6 ConclusionIn this chapter, we have proposed a chemo-mechanical model for neu-trophil activation in response to mechanical deformation, and simulated thefluidization, polarization and activation of a neutrophil as it traverses anarrow microfluidic channel. The model integrates insights and hypothesesfrom previous experiments into a coherent theoretical framework, and pre-dicts the salient features of how a neutrophil senses mechanical stimulationand responds by remodelling its cortex. Thus, the process constitutes a casestudy on modelling cell dynamics based on the dual pillars of mechanics andchemical kinetics. We may summarize the main results as follows:• Fluidization. As the neutrophil deforms and enters the narrow channel,its cytoskeleton starts to melt at the cell front where it sustains severemechancial deformation. The melting front propagates toward the rearof the cell, and by the time the bulk of the cell is inside the channel,the entire cytoskeleton has fluidized.• Polarization. Melting of the cytoskeleton eliminates the cortex tensionthat inhibits Rac activity on the membrane. Thus, inactive Rac in thecytosol becomes active Rac on the membrane. If the kinetic rate of Racturnover is high relative to the reciprocal of the entry time, a polarizedRac distribution prevails, and persists by a wave-pinning mechanismeven after the cell exits the confining channel.• Activation. The polarized Rac distribution induces cortex recovery inthe front as well as a protrusion force that generates a pseudopod.Thus the neutrophil becomes activated.Thus, the model captures the main features of the experiment from a setof reasonable assumptions about the chemical kinetics, signalling pathways,and hydrodynamics. In particular, the results have confirmed the hypothesis[237, 238] that neutrophil activation depends on the competition of two timescales. If the kinetic time is shorter than the cell entry time, Rac has time714.6. Conclusionto polarize and the cell is activated. If the kinetic time is too long, themechanical deformation is relieved before polarization and activation canoccur, and the cell simply returns to its spherical inactive state. Based onthis competition, we have constructed a phase diagram from numerical data,and the critical condition for cell activation is comparable to experimentalobservations. The shape of the activated cell and the polarized distributionsof Rac, cortical modulus and protrusion force also agree qualitatively withexperiments.Of necessity the model introduces various simplifications, and these mayhave limited the model’s ability to capture more refined features of the realprocess. For example, the predicted “pseudopod” consists of the entire frontportion of the cell, whereas in reality, pseudopods take on smaller dimen-sions and varied shapes and locations. The model simplifications fall intotwo categories, mechanical and chemical. Mechanically we model the cellin 2D, and represent the plasma membrane and the underlying cortex us-ing a single mechanical element: linearly elastic springs connected at nodesalong the cell’s outline. As such, we neglect the details of cortex break-age and destruction, and indicate the cortical integrity only by the corticalmodulus. Finally, the cytosol is assigned the same viscosity as the outsidefluid. This simplifies the immersed-boundary formalism and facilitates thecomputation. In treating the chemical kinetics of polarization, we have ne-glected the intricate relationship among Rho GTPases, and instead used asingle protein, Rac, to describe the polarization process. Because of this,the inhibition of membrane Rac by cortical tension, the promotion of corti-cal reinforcement by Rac, and the development of a protrusion force are allmodeled ad hoc.These assumptions and simplifications can be viewed as motivations forfuture work. Validation of the model assumptions calls for future experi-ments that can measure protein and force distributions with spatial resolu-tion. More refined models can systematically remove the limitations, anddetermine the extent to which they have distorted the true dynamics of theprocess.72Chapter 5Fluidization of cells underdynamic strain5.1 IntroductionConventional wisdom, informed by polymer physics, holds that stretch-ing causes cells to stiffen. Known as the reinforcement effect, the stiffeningis manifested by increase of one or more of the following experimentallymeasured properties: the elastic modulus of the cytoskeleton [213], cell stiff-ness probed by micro-bead displacement [133], traction force on deformablesubstrate [59], and formation of new focal adhesion sites [60]. More recentexperiments by Fredberg and coworkers [18, 108, 109, 213], on the otherhand, demonstrated an apparently opposite response known as fluidization,indicated by a drastic decrease of the elastic modulus and disintegration ofthe cytoskeleton of the cells.To rationalize this apparent contradiction, several studies have compareddifferent modes of imposing mechanical deformation. In a “stretch-hold”(SH) maneuver, the cell is stretched to a prescribed strain and then heldat the deformed state. This mode of deformation engenders no controversy;all studies reported cell stiffening similar to a polymer network under exter-nal stress [149]. In a stretch-compress (SC) maneuver, the cell is stretchedto a prescribed strain ε0 and then compressed so as to yield a net strainof zero. Following such a maneuver, the cell stiffness is seen to drop by aconsiderable amount, about 50% for ε0 ∼ 0.1 [18, 109]. This is fluidization.Afterwards, the cell gradually recovers its equilibrium stiffness over some 5minutes. Gavara et al. [59] and Chen et al. [18] further demonstrated thatthe fluidization happens during the compression stage of the maneuver; theinitial stretch in fact increases cell stiffness as in the SH maneuver. Remark-ably, a compress-stretch (CS) maneuver, where a cell is first compressed andthen returned to zero strain, does not elicit fluidization [18].Taken together, these studies have established a history-dependence ofthe mechanical response of the cell. But what is its mechanism? Is it735.1. Introductionpurely mechanical, akin to the thixotropy of complex fluids? Or is it uniqueto biological cells, arising from biochemical signals that direct structuralremodeling? What causes the stark “asymmetry” between the CS and SCmaneuvers in fluidizing the cell?One expects the mechanical response, be it stiffening or fluidization, to bedirectly related to structural changes in the cytoskeleton. The most promi-nent entity in these adherent cells is the stress fibers, which are bundles ofactin filaments cross-linked by myosin, α-actinin and other cross-linkers, fre-quently emanating from focal adhesions. Earlier experiments [23, 182, 183]have shown that compression of a large enough strain causes stress fibersto buckle and disintegrate, in accord with the fragility of actin filamentsunder compression [200]. This suggests that the disassembly of the stressfibers causes fluidization. In support of this idea, Chen et al. [18] reporteddisassembly of the cytoskeleton during the compression phase of the SC ma-neuver. Besides, cyclic stretch-unstretch deformation of cells causes disas-sembly of stress fibers oriented along the direction of deformation and theirprevalence in the orthogonal direction [17, 96, 108, 212]. However, howdoes this mechanism discriminate between SC and CS maneuvers? Why isthe initial stretch prerequisite to subsequent disassembly of the stress fibersduring compression?To our knowledge, two models have appeared on cell fluidization, bothbased on the dynamics of semi-flexible polymer networks. First, Morozovand Pismen [149, 150] showed that mechanical stretching may cause ananisotropic redistribution of myosin in the actin network, thereby lower-ing its modulus and causing fluidization. However, this argument based onanisotropic prestress does not distinguish between the SC and CS maneu-vers, and thus does not speak to the dependence of fluidization on defor-mation history. Besides, the mechanism of disassembly for actin filamentsis rupture under high tension. This differs from the conventional view ofbuckling under compression [18, 23, 200]. Second, Wolff et al. [233] built awormlike-chain model to explain the rheology of reconstituted actomyosinnetworks under large-strain oscillatory shear. The network is predicted tostiffen after the first cycle of strain, and then gradually soften with subse-quent cycles due to breaking of transient bonds between polymer chains.This bears some resemblance to cell fluidization. But the latter happensafter a single stretch-compress cycle [18, 108]. There is also no discussion ofthe contrast between SC and CS maneuvers.Models of cyclic stretching of cells [17, 96, 108, 212] are also potentiallyrelevant to cell fluidization. These typically attribute stress fiber disassemblyto the detachment of myosin motors, with a probability based on the sliding745.2. Model for stress fibersvelocity of the motors on the actin cable and the tensile force in it [212].This predicts disassembly of the stress fibers during the stretching phase ofthe cycle, as opposed to the compression phase [18, 59]. In addition, Chenet al. [18] suggested an analogy between cell fluidization and the behavior of“catch bonds” [17, 129], whose life time increases under stretching. Althoughcommon in focal adhesions, catch bonds have never been reported in thecytoskeleton [18], and its molecular basis is unclear at present.The objective of this chapter is to explain the dependence of cell fluidiza-tion on the stretching protocol, more specifically the contrast between thestretch-compress (SC) and compress-stretch (CS) maneuvers. We couple asimple mechanical model for the stress fiber with the crossbridge cycle ofthe myosin motors, and demonstrate that in the SC maneuver, the initialstretching produces a condition in which a subsequent compression causes alarge fraction of the myosin to detach from the actin filaments. This is takenas the precursor to the disassembly of the stress fiber and the cytoskeleton,i.e. fluidization. The CS maneuver, on the other hand, does not producethe same condition.5.2 Model for stress fibersDespite the various puzzles mentioned above, one fact seems well-establishedby experiments: the fluidization corresponds to the disassembly of the cy-toskeleton [18, 108, 109, 213]. This is the physical entity underlying allpreviously used measures of fluidization—losses in cell modulus, tractionforces on the substrate and focal adhesions. Following the spirit of numer-ous prior studies [17, 96, 212], we will focus on the dynamics of a single stressfiber (SF). Such a dramatic simplification, down from the complexities of theactomyosin network or even the whole cell, allows us to capture the essen-tial features of the fluidization process in a minimalist fashion. Such anSF model should capture three key effects. First, the SF is a viscoelasticsolid that sustains a finite prestress in equilibrium [98, 149]. Second, theSF has a myosin contractile element that imposes the prestress in equilib-rium. During mechanical deformation, the myosin apparatus should providea relaxation mechanism that adjusts the “unloaded reference length” of thestress fiber [98]. These two considerations have led us to the Kelvin-Voigt-Myosin (KVM) model depicted in Fig. 5.1, with the myosin unit in serieswith the spring-dashpot assembly. The third key effect is the power-strokeof the myosin crossbridge, which will be discussed shortly.The KVM model is based on previous models for the SF [8, 17, 89, 98,755.2. Model for stress fibersFigure 5.1: Our model for the stress fiber has a Kelvin-Voigt element con-nected in series to a myosin apparatus. The two ends of the stress fiberare attached to the substrate through focal adhesions, which are assumedpermanent in our model.178] and smooth muscle cells and strips [7, 14]. Strictly speaking, it onlyrepresents one sarcomeric unit of the SF. As is usually done [17], we tacitlyassume that all sarcomere units are activated simultaneously. Besser andSchwarz [8] used a Kelvin-Voigt element in parallel with a myosin appara-tus to study the dynamic coupling between the stress fiber and the focaladhesion. Laser ablation of stress fibers, however, argues against a parallelconfiguration in favor of one with the elastic and myosin elements in series[178]. Later, the series configuration was used to study the dynamics of SFunder cyclic stretch [89, 98]. To reflect the viscoelasticity of the SF [110], wehave further added an viscous element to this model to arrive at the KVMconfiguration of Fig. 5.1. A similar model has been used by Chen et al. [17]to study SF reorientation under cyclic stretch. Thus, the total length of thestress fiber is L = LKV + Lm, with LKV and Lm being the lengths of theKelvin-Voigt element and the myosin element, respectively. The myosin andKelvin-Voigt elements share the same tension, which is also the tension ofthe SF as a whole:σ = G(LKV − L0) + ηL˙KV , (5.1)where G is the elastic modulus and η is the viscosity of the Kelvin-Voigtelement, and L0 is the resting length of the Kelvin-Voigt element. The dotindicates time derivative.Myosin plays a key role in our model, partly because its detachment enmasse from the actin filament will be taken to mark the onset of disassemblyof the stress fiber, which leads to fluidization of the cell. The experimentalevidence for such a criterion for SF disassembly has been reviewed by Kau-nas and Deguchi [96]. In particular, Trepat et al. [213] and Chen et al. [18]765.2. Model for stress fibersFigure 5.2: (a) The swinging crossbridge model adapted from Spudlich [201].(b) A simplified three-stage model represents the interaction between themyosin and actin filaments, with m0, m1 and m2 denoting myosin in thedetached, attached pre-stroke and attached post-stroke states. These statescorrespond roughly to those marked in (a).have ascribed the disruption of the cytoskeleton to the breakage of cross-linkers, especially the detachment of the myosin crossbridge. Gavara et al.[59] attributed cell fluidization after the SC maneuver to “a fall in active con-tractile tension caused by stretch-induced actomyosin detachment.” Matsuiet al. [131, 132] showed that high levels of MgATP reduce the number ofactomyosin crossbridges, and promote rapid disassembly of the SF. Further-more, there exists a clear connection between the loss of tension in the SFand its disassembly [164]. Compression of SFs reduces the tensile prestressin them, and causes the SF to disassemble, or if the strain rate is sufficientlyhigh, to buckle [23, 182, 183, 200]. These ideas inform the third componentof our model, the cycle of the myosin crossbridge. Our goal is to show thatthe SC maneuver generates a loss of tension in the KVM model of the stressfiber, such that myosin motors detach from the actin filaments and causedisassembly of the stress fiber. A CS maneuver, on the other hand, will nothave the same effect.To capture the kinetics of the myosin motors, we start from the well-known swinging crossbridge model (Fig. 5.2a) [201]. First, myosin carryingan ADP and a phosphate binds to an actin filament in the configurationmarked with m1. Release of the phosphate prompts the light-chain-bindingregion of myosin to swing through an angle, providing a working strokeof 5–15 nm. This motion, known as the “power-stroke”, is where the me-chanical work of contraction is done. Next, the ADP is released, and thebinding of a fresh ATP causes the myosin to detach from the actin filament.Hydrolysis of the ATP then returns the myosin to the “cocked” conforma-775.2. Model for stress fiberstion and primes it for the next cycle. For our purpose, we represent theswinging crossbridge cycle by three steps connecting three distinct myosinstates: the detached, the attached pre-stroke and the attached post-strokestates (Fig. 5.2b). This can be viewed as a simplified version of the pop-ular crossbridge model for smooth muscle cells [19, 41, 69, 135, 227], withthe unphosphorylated myosin motors disregarded. The amount of myosinin each of the states are, respectively, m0, m1 and m2. Neglecting reversereactions, we can describe the myosin kinetics by the following equations:mt = m0 +m1 +m2, (5.2)m˙1 = k01m0 − k12m1, (5.3)m˙2 = k12m1 − k20m2, (5.4)where k01, k12 and k20 are the rate constants. We assume that k01 and k12remain constant during the crossbridge cycle. The rate of myosin detach-ment k20, on the other hand, is known to decrease with mechanical load onthe motor [107, 220, 221]. In particular, a myosin motor may “stall” on anactin filament when the tensile force σ reaches a critical stalling force σ0.Following Veigel et al. [220, 221], we assume an exponential form for k20:k20 = Koffe−k(σ−σ0). (5.5)This relationship will prove to be important to our model, as it representsthe sensitivity of the myosin duty ratio to external forcing on the stress fiber.Similar to the models of Kaunas et al. [96, 98], our myosin apparatusactively regulates the length of the stress fiber. In equilibrium, the myosinmotors walk until the elastic tension in the SF equals its stalling force σ0.Further stretching or compression will cause the motors to slide or walk inan attempt to return the SF to “tensional homeostasis” [15, 98]. The re-lationship between the contractile speed of the myosin apparatus and themechanical load is commonly represented by Hill’s law [7, 78, 191]. In SFmodels, however, a linearized version is often used to avoid algebraic com-plexity [17, 21, 98, 178, 202, 231]: v/Vm = (1 − σ/σ0), where Vm is theupper bound of myosin speed attained under zero load. Furthermore, theexternal load not only changes the myosin speed directly as above, but alsoindirectly by modifying the myosin duty ratio, and therefore the fraction ofmyosin motors that are actively “walking” and contributing to the velocityv. To account for this effect, we modify the velocity-load relationship asv = −L˙m =m1mtVm(1−σσ0). (5.6)785.3. Model parametersThis modification is based on two considerations. First, during the cross-bridge cycle [201], only attached pre-stroke motors (m1) are capable of pro-ducing contraction. After the power-stroke, the myosin molecules are stuckin the rigor state pending arrival of the next ATP, which will cause themto detach. Therefore, it seems reasonable to make v proportional to m1.Second, Warshaw et al. [230] have designed an in vitro motility assay inwhich a fraction of the myosin motors have an unphosphorylated light chainand hence remain inactive. They found that the unphosphorylated myosinmotors do not generate force but act as a load to slow down the contrac-tion that is driven by the cycling phosphorylated myosin motors. Drawingan analogy with our model, we consider the motors not actively walking,i.e., m2 and m0, to act as dampers to the contraction. Hence we haveput the total myosin number mt in the denominator. Warshaw et al. [230]suggested a sigmoidal dependence of the sliding speed on the fraction ofphosphorylated myosin. We adopt a linear relationship for simplicity. Notethat our model assumes the total number of myosin on the stress fiber mtto be fixed. Micropipette aspiration experiments [122, 123, 198] suggestthat external forcing may recruit myosin motors from the cytoplasm to theactomyosin cortex. Our simple model disregards such spatial transport.Equations (5.1–5.6) form a complete description of our model stress fiber.When a prescribed deformation is imposed on the total length L, theseequations can be time-stepped for the evolution of the tension σ, the myosinpopulation in each of the three states m0, m1 and m2, and the lengths ofthe KV and myosin elements. One feature of this model might be counter-intuitive. The instantaneous tension σ(t) is not directly given by the numberof attached myosin molecules. The latter controls the sliding or walkingvelocity of the myosin apparatus, which then modulates the relaxation ofthe KV element. It is through these two intermediaries that myosin controlsthe tension σ.5.3 Model parametersTo introduce the parameters of the problem, it is convenient to considerfirst the equilibrium or homeostatic state of the KVM apparatus. Givenany total length L for the apparatus, the myosin element will contract (bysub-stall walking) or expand (by super-stall sliding) to adjust its length andthe length of the KV element such that in the equilibrium state, the tensionin the apparatus is exactly the stalling force σ0. This determines the “pre-strain” of the KV element to be εp = σ0/(GL0). In equilibrium, the myosin795.3. Model parametersin different states can be easily calculated from Eqs. (5.2–5.4):m0 =k12Koff(k01 + k12)Koff + k01k12mt, (5.7)m1 =k01Koff(k01 + k12)Koff + k01k12mt, (5.8)m2 =k01k12(k01 + k12)Koff + k01k12mt. (5.9)Note that the equilibrium myosin length Lm and total length L are notunique but depend on the total strain imposed on the KVM apparatus. AsLm enters the dynamics only through L˙m (cf. Eq. 5.6), we set Lm = 0 atthe equilibrium state without losing generality. For this reason, the strainfor the subsequent deformation is given in terms of the resting length L0instead of L.Starting from equilibrium, the SC maneuver is imposed by stretching theKVM model at a constant speed to a prescribed total extension of L0ε0, andthen compressed at the same constant speed to the original length. Denotingthe duration of SC maneuver by T , the maneuver can be represented by thestretching velocity:L˙ ={ 2L0ε0T if t < T/2−2L0ε0T if t > T/2(5.10)This amounts to a triangular wave form of the strain. The CS maneuver isthe reverse of the above.Table 5.1 lists all the parameters of the problem, their values as adoptedhere and their sources. The determination of several of the parametersrequires some explanation. We define the duty ratio γ as the fraction of timethat a myosin motor spends being attached to the actin filament during acrossbridge cycle in equilibrium [22, 177]. Thus, Eqs. (5.8) and (5.9) give usγ =k01k12 + k01Koff(k01 + k12)Koff + k01k12. (5.11)In Table 5.1, we have adopted the duty ratio γ = 0.8 appropriate for themyosin IIB isoform. For the IIA isoform, the duty ratio is much lower,around γ = 0.1 [22]. Myosin IIB binds to actin filaments for much longertime under resistant stress [107], and is thus more effective in maintainingtension in the stress fibers and more sensitive to mechanical strains [226].While both myosin IIA and IIB contribute to stress fibers in adherent non-migrating cells [105], they assume spatially polarized patterns in migrating805.3. Model parametersTable 5.1: Parameters for the KVM stress-fiber model.Parameter meaning Value and referencesγ myosin duty ratio 0.8 [22, 177]k01 myosin attachment rate 4 (s−1)k12 phosphate release rate 1.5 (s−1) [177, 226]Koff myosin detachment rate 3 (s−1) [107, 177, 220, 226]k mechanical sensitivity (Eq. 5.5) 0.6 (nN−1) [31, 32, 107, 229]L0 resting length of KV unit 1 (µm) [180]G elastic modulus 50 (nN/µm) [8, 17, 32, 110]η viscosity 200 (nN·s/µm) [17, 110]εp pre-strain 0.1 [23, 121]σ0 myosin stall force 5 (nN) [32]Vm characteristic contractile velocity 0.075 (µm/s) [17, 94, 95, 98, 132]cells. IIB is concentrated at the rear of the cell in the contractile stressfibers, whereas IIA is predominantly found in the protrusive actin structuresat the leading edge [106, 223]. This suggests distinct functions of IIA andIIB in forming different actomyosin structures. For its prominent role inmaintaining the contractile stress fibers at the rear of migrating cells [222,223], IIB appears more relevant to the stress fibers considered here, andhence its duty ratio is used.Using measured values of k12, Koff and γ, Eq. (5.11) allows us to deter-mine the myosin attachment rate to be k01 = 4 s−1. Thus, in the equilibriumstate, the myosin motors are distributed in the three states as m0 = 0.2mt,m1 = 0.533mt and m2 = 0.267mt. The mechanical sensitivity parameterk is determined from the measurement of Kova´cs et al. [107] on a singlemyosin head and an estimation of the number of myosin motors in a stressfiber [31, 32].The characteristic contractile velocity Vm is estimated from the contrac-tion rate of isolated, load-free SFs [94, 95, 132]. In our model, the distribu-tion of myosin motors under zero load is such that m1/mt = 0.714. ThroughEq. (5.6), the measured SF contraction rate of 0.01–0.2 s−1 corresponds to0.014 ≤ Vm ≤ 0.28 µm/s. Vm = 0.075 µm/s is chosen as a representativevalue. As discussed later, this choice is also based partly on fitting the ex-perimental data in Fig. 5.7 (see also Fig. 5.8). Thus, Vm is the sole fittingparameter of the model.The ODE system is solved by using a fourth-order Runga-Kutta schemein MATLAB. Each simulation starts from the equilibrium state with a pre-815.4. ResultsFigure 5.3: Tension in the stress fiber under the SC (a) and CS (b) maneu-vers.strain of εp = 0.1. The time step is set to be ∆t = 0.004 s. Numerical testsshow that this ensures sufficient temporal resolution; refining ∆t furthermorebrings negligible changes to the results.5.4 Results5.4.1 Responses to the SC and CS maneuversTo study cell responses to the SC and CS maneuvers, Chen et al. [18]applied a strain magnitude of ε0 = 0.1 over T = 4 s, with a strain rateof ε˙0 = 0.05 s−1. The SC maneuver induces cell fluidization, while the CSmaneuver does not. In our simulation, we similarly apply the two maneuversat ε0 = 0.1 and ε˙0 = 0.05 s−1. The rest of the parameters are as listed inTable 4.1.For the SC maneuver (Fig. 5.3a), the tension σ jumps up at the start,as the finite strain rate instantaneously elicits a finite force from the viscouselement. Then σ continues to increase thanks to the stretching of the elasticspring. The rate σ˙ itself increases in time. This is because the stretchreduces the number of myosin motors m1 and thus the sliding speed of themyosin element (cf. Eq. 5.6). Consequently, the KV element is forced toincrease its rate of stretching. Hence the faster rise of σ. As the stretchingis reversed at t/T = 0.5, the tension σ abruptly drops due to the reversal inthe viscous force, to a negative value in this case. Then σ recovers in timeto a positive value that is roughly 40% its initial equilibrium value.Upon starting the CS maneuver (Fig. 5.3b), the tension σ drops instanta-825.4. Resultsneously due to sudden imposition of a negative strain rate. It then maintainsa roughly constant level through the whole compression phase. In this pro-cess, the myosin motors walk at a more or less constant speed, and theKV element sustains a roughly constant rate of compression. The tensionσ reflects mostly the viscous contribution, similar to the compression stageof the SC maneuver. At the start of the subsequent stretching phase, thetension suddenly jumps up and then gradually increases due to elongationof the Kelvin-Voigt element. For the parameters used, our model predictsthat the SC maneuver produces a negative (compressive) force within thestress fiber while the CS maneuver does not.Now we can examine the evolution of the three myosin populations dur-ing the SC and CS maneuvers. For the SC maneuver (Fig. 5.4), the risingtension during the stretch phase lowers the myosin off-rate k20 (cf. Eq. 5.5)and suppresses myosin detachment from actin filaments. This causes ac-cumulation of myosin in the m2 stage and depletes the m1 population(Fig. 5.4a). But the total number of attached myosin ma = m1 + m2increases moderately during the stretch phase (Fig. 5.4b). At the end ofstretch (t = 0.5T ), m2 myosin makes up about 90% of all the myosin mo-tors. More importantly, the low level of m1 means that now the myosin ele-ment has a very slow walking speed. Upon reversal of the deformation, themyosin element cannot contract quickly enough to absorb the compressionimposed on the entire stress fiber. Thus, the negative strain rate is mostlysustained by the KV element, which produces the negative σ of Fig. 5.3(a).The compressive force significantly increases the off-rate k20 and results ina rapid detachment of m2 myosin. The total number of attached myosinma = m1 + m2 undergoes a sharp drop as well (Fig. 5.4b). This scenariocorresponds to the “loss of tension” that is the immediate trigger for stress-fiber disassembly in prior experiments and models [23, 164, 182]. In ourmodel, the newly detached m0 motors reattach at the fixed rate k01, andthus ma recovers to some extent during the compression (Fig. 5.4b). In real-ity, the cytoskeleton rapidly disintegrates, and only reforms several minutesafter the SC maneuver has been completed [18].The reaction of the myosin populations to the CS maneuver is entirelydifferent. Upon the start of compression, m2 rapidly decreases while m1increases moderately (Fig. 5.5a). The increasing m1 enhances the walkingspeed of the attached myosin and the relaxation rate of the myosin appa-ratus, which helps to alleviate the compression being sustained by the KVelement. Consequently, the drop in the tension σ is much smaller than in thecompression part of the SC maneuver (Fig. 5.3). In particular, it remainstensile throughout the compression phase. As a result, the compression only835.4. ResultsFigure 5.4: Evolution of myosin populations under the SC maneuver. (a)m1 and m2; (b) the total attached myosin ma = m1 +m2.causes a moderate loss of m2 and an even smaller drop in the total attachedmyosin ma (Fig. 5.5b). The subsequent stretching recruits more myosinonto the actin filaments, elevating m2 and ma and raising the tensile forceσ. Note that the contrast between SC and CS maneuvers illustrated abovedepends on the high duty ratio of myosin IIB. When the low duty ratio ofIIA is used, the SC-CS difference is much reduced, and the force in the stressfiber stays tensile in both cases.At the end of the SC maneuver (t = T ), the force σ recovers instanta-neously to a value slightly above the stalling force σ0 (not shown in Fig. 5.3).This can be rationalized as follows. The compression phase enjoys a larger845.4. ResultsFigure 5.5: Evolution of myosin populations under the CS maneuver: (a)m1 and m2; (b) the total attached myosin ma = m1 +m2.m1 population than the stretching phase (Fig. 5.4a), and thus faster myosinwalking (Eq. 5.6). As a result, the elastic element has compressed less thanit has expanded, ending up in a stretched state (LKV > L0) that produces atension above σ0. This relaxes toward the equilibrium state in time, fallingwithin 1% of σ0 in 1.5T . Similarly, the end of the CS maneuver sees anabrupt drop of σ to a level above σ0, followed by a gradual relaxation to-ward equilibrium.In summary, our model accounts for the distinct reaction of the stressfiber to SC and CS maneuvers from the dynamics of the different myosinpopulations. Essentially, the stretching in SC causes a great number of855.4. Resultsmyosin motors that have completed the power stroke to stay attached ontoactin filaments. Furthermore, this myosin distribution reduces the speed atwhich the myosin element can relax. Both factors conspire to produce therapid detachment of a large fraction of the myosin motors when the compres-sion phase starts. The consequence is disassembly of the stress fiber. Thus,the stretching establishes a precondition for the subsequent detachment ofmyosin upon compression. For the CS maneuver, such a precondition isabsent. As a result, the initial compression produces a much milder loss ofattached myosin and no fluidization.The SC-CS asymmetry discussed above is quite robust and insensitiveto certain model assumptions. For example, it still arises if the Kelvin-Voigt element is replaced by a Maxwell element, as long as the stretch putsa large enough force on the myosin apparatus that markedly reduces themyosin sliding velocity. Similarly, if the nonlinear Hill’s law is used in placeof the linear velocity-force relationship (Eq. 5.6), a slower sliding velocityprevails under super-stall forcing. Thus, during the stretching phase, themyosin apparatus will sustain a larger force. This accumulates a largernumber of attached myosin motors, which will subsequently detach en masseupon start of the compression phase. Hence, the Hill model will accentuatethe fluidization response. Finally, the SC-CS asymmetry can be producedas long as the myosin velocity increases with the proportion of attachedmyosin m1/mt. Thus, the myosin apparatus becomes more rigid during thestretching phase, with the depletion of m1, setting up the condition for theloss of tension upon reversal of the strain rate. The linear dependence onm1/mt in Eq. (5.6) is not essential but adopted for its simplicity.5.4.2 Stretch-hold (SH) maneuverThe stiffening of the cell in response to an SH maneuver has been demon-strated in several experiments [59, 139, 167, 213]. All experiments havereported the same general features. The tensile force σ rises during thestretching. Upon cessation of stretching, σ experiences a prompt decreasefollowed by a more gradual relaxation toward σ0. Individual studies differin quantitative details. For instance, Pourati et al. [167] applied strains of2.5% and 5% in a short time period (< 10 s) on adherent endothelial cells,and recorded a 15% and 30% increase in cell stiffness within 70 seconds ofthe cessation of the stretch. Mizutani et al. [139] have reported a roughly40% increase in Young’s modulus for fibroblasts after an 8% stretch. Thenthe modulus decreases gradually in the following 2 hours. Trepat et al. [213]applied a 10% stretch on human airway smooth muscle cells over 2 seconds,865.4. ResultsFigure 5.6: Evolution of the tension σ during an SH maneuver modeled afterthe experiment of Trepat et al. [213]: a 10% strain is applied over T/2 = 2s and then maintained indefinitely. The characteristic time T = 4 s. ∆σ isthe excess tension after cessation of the stretch, and td is the time for ∆σto relax by 50%.and reported a 50% increase in their stiffness immediately afterwards. Inthe ensuing “hold”, it took more than 10 min for the cell to relax to its orig-inal stiffness. Gavara et al. [59] observed apparently stronger effects thanthe other groups, documenting an increase of 56% and 77% in cell-substratetraction after a 5.5% and 11% stretch, respectively. The increased tractionwas measured within 2 minutes of holding the stretch, but they did not holdthe stretch longer to explore the long-term relaxation toward equilibrium.Inspired by these experiments, we have carried out stretch-hold simula-tions using our KVM model. The stretch is implemented at a constant ratefor a time T/2 (cf. Eq. 5.10), and then the stretched state is maintainedindefinitely. Figure 5.6 shows the model prediction of the tension σ(t) us-ing the same parameters as in the SC maneuver (Fig. 5.3a). As expected,reaction to the initial stretch is identical to that in the SC maneuver, withthe tension increase in time to approximately 3.2σ0 at t/T = 0.5. Afterthe stretch stops, the tension abruptly drops to about 1.5σ0 due to loss of875.4. ResultsFigure 5.7: Comparison of the excess tension, upon cessation of stretch,between the model prediction and experimental data.the viscous force. Afterwards, σ gradually relaxes toward the homeostaticlevel σ0. The viscoelastic response observed here is comparable to that offibroblast in uniaxial stretching [209], where a similar Kelvin model is usedto interpret the data. Two quantities of interest are ∆σ, the excess tensionabove the homeostatic level immediately after the cessation of stretch, andthe “half-life” td for the subsequent relaxation, the time for the excess ten-sion to drop by 50%. These will be compared with experimental values inthe following.Figure 5.7 compares ∆σ as a function of the strain magnitude with ex-perimental data [59, 139, 167, 213]. Two caveats must be noted. First,only Trepat et al. [213] managed to capture the cell stiffness more or lessimmediately upon cessation of stretching, thanks to their optical magnetictwisting cytometry technique. The other experiments recorded such dataa short time, roughly 1 to 2 minutes, afterwards. Thus, only Trepat et al.[213] is strictly comparable with our model prediction. Second, most of theexperiments measured cell modulus by microrheometry and magnetic twist-ing cytometry. Since the elastic modulus of biological cells is proportionalto the contractile force in the actomyosin network [228], the percentage ex-885.4. Resultscessive modulus is comparable with the percentage excess tension. Withthese in mind, Fig. 5.7 shows that all the experimental data lie close to themodel prediction using Vm = 0.075 µm/s. In particular, there is quantitativeagreement with the result of Trepat et al. [213]. Our model also capturesthe rising trend in the data; a larger strain ε0 causes the elastic element toextend more at the end of stretching.Unfortunately, the model fails to capture the time scale of the subse-quent relaxation of the excess tension. The calculation in Fig. 5.6 matchesthe experimental conditions of Trepat et al. [213] in that a 10% strain is ap-plied linearly over 2 seconds. The relaxation has a half-life of td = 1.9 s, andas a whole lasts about 16 seconds. Experimentally, on the other hand, td isaround 400 s and the whole relaxation lasts about 17 minutes. Thus, therelaxation is too fast in the simulation by two orders of magnitude. One pos-sible explanation is that Trepat et al. [213] measured whole-cell relaxationwhile the model deals with a single SF. Thus, slow diffusive processes ofre-equilibrating a large number of SFs are missing from the model. Besides,the model uses a linear force-velocity relation (Eq. 5.6) instead of the morerealistic hyperbolic Hill’s law [31, 78]. Under super-stall force, the formerpredicts a much larger sliding velocity than the latter. As the SH maneuverresides entirely in the super-stall regime, the relaxation would proceed muchmore rapidly than in reality.Aside from the SH maneuver, several experiments have explored thereaction of adherent cells to the stretch-compress [18, 59, 109], compress-hold[23, 182, 183] and cyclic stretch [89, 96–99, 212] maneuvers. The most salientfeature of these experiments is the critical conditions for the buckling anddisassembly of the stress fibers. It is possible to compare model predictionsquantitatively with the measured critical strains and strain rates, and wereport such comparisons next.5.4.3 Critical conditionsPrior experiments have subjected adherent cells to different deforma-tions, and recorded the critical strains and strain rates that are requiredfor the disassembly of stress fibers and fluidization of the cells. Thus itis interesting to explore the model predictions of these critical conditions,and compare them with experimental measurements. In our Kelvin-Voigt-Myosin (KVM) model, we will take the complete loss of tensile force in theSF, i.e. σ = 0, to be the critical condition. Such a choice is motivated byexperimental observations that a stress fiber buckles and disintegrates afterits pre-tension is relaxed by compression [23, 164, 182]. It is also consistent895.4. Resultswith two microscopic effects. The first is the highly asymmetric load re-sponse of actin filaments: they can sustain considerable tension but buckleeasily under pico-newtons of compression [200]. The second is the enhanceddetachment rate of myosin motors due to loss of tensile force [107, 220] (cf.Eq. 5.5).Motivated by experimental data, we have chosen to construct a phasediagram in the strain–strain rate plane. Using the stretch-compress (SC)maneuver of Eq. (5.10), we carry out a series of simulations by varying eitherthe total strain ε0 or the strain rate ε˙0 while holding the other constant. Asε˙0 = 2ε0/T , this is accomplished by adjusting the period T . The outcome isrecorded as either fluidization or non-fluidization, depending on the sign ofthe minimum tension. The boundary between the two regimes is depictedin Fig. 5.8(a). Experimental observations suggest that fluidization occurswhen either ε0 or ε˙0 is high enough. Thus, one may have anticipated a sortof compensation between them, and a negative slope for the boundary inthe (ε0, ε˙0) plane. The model shows, however, a non-monotonic boundarywith a positive slope for large strain.The key to understanding this behavior is to recognize that in our KVMmodel, the force σ in the SF is affected by the relaxation of the myosinapparatus. For a large strain ε0, σ grows so much during the stretching phase(cf. Fig. 5.3a) that at its end (t = T/2), almost all the m1 myosin motorshave been converted to m2 (cf. Fig. 5.4a). As a result, myosin walkingand relaxation become negligible according to Eq. (5.6), and the myosinapparatus is essentially rigid. The subsequent compression, therefore, issustained almost entirely by the Kelvin-Voigt (KV) element. In this regime,a larger ε0 implies a larger elastic tension from the KV element at the end ofstretch, which will protect the SF from disassembly during the subsequentcompression. Thus, a faster strain rate ε˙0 is required for disassembly atlarger strain ε0, and the slope of the boundary should be positive. To bemore precise, we estimate the minimum tension at the start of compressionas σmin = GL0(εp+ε0)−ηL0ε˙0, where we have neglected the strain sustainedby the myosin apparatus. Setting σmin to zero then gives us a criticalcondition for fluidization that is a straight line of positive slope G/η inthe (ε0, ε˙0) plane. This explains the right portion of Fig. 5.8(a); the actualslope of the boundary is 0.173 s−1, reasonably close to the expected G/η =0.25 s−1. The difference is due to the sliding of the myosin motors. Torationalize the negative slope for small ε0, we recognize that with decreasingε0, the myosin apparatus is more flexible with faster relaxation at the startof compression. Thus, the KV element will suffer less strain, and fluidizationis achieved only with greater strain rate ε˙0. Hence the boundary takes on905.4. ResultsFigure 5.8: Critical conditions for SF disassembly depicted on the (ε0, ε˙0)plane. For the experimental data, open symbols correspond to fluidization,filled ones non-fluidization, and half-filled ones the critical condition. (a)The SC maneuver; (b) The CS and CH maneuvers. Model prediction at alower Vm = 0.025 µm/s is also shown for comparison.a negative slope at the left portion of Fig. 5.8(a). Incidentally, most priorexperiments fall within this portion (ε˙ < 0.1), and thus exhibit fluidizationwhen either ε0 or ε˙0 is high enough.The critical conditions for compress-stretch (CS) and compress-hold (CH)maneuvers also present intriguing features (Fig. 5.8b). As the minimum ten-sion (or largest compressive force) occurs during the compression phase ofboth CS and CH, the critical condition is the same for both maneuvers.Upon onset of compression, σ drops instantaneously to a local minimumthanks to the viscous force (cf. Fig. 5.3b). As σ falls below σ0, the myosinapparatus starts to contract to alleviate the compression on the KV ele-ment, at a speed of (m1/mt)Vm(1 − σ/σ0) = 0.533Vm(1 − σ/σ0) accordingto Eq. (5.6). Note that the equilibrium myosin ratio m1/mt = 0.533 holdsas there has not been enough time for the myosin populations to evolve.Thus, for a small strain ε0, the critical condition for fluidization amounts toσmin = GL0εp − η(L0ε˙0 − 0.533Vm) = 0. We have ignored the elastic forcedue to compression because the strain ε0 is small. This predicts a criticalcondition for small strain ε0:ε˙0 = 0.533VmL0+Gηεp. (5.12)For a large strain, on the other hand, the critical condition of σ = 0 amountsto the myosin apparatus contracting at the zero-load speed of 0.714Vm (note915.4. Resultsm1/mt = 0.714 at zero load), with the KV element completely relaxed. Thisgives a large-strain asymptote for the critical condition: ε˙0 = 0.714Vm/L0.Both limits are on the strain rate alone, regardless of the strain. Thus, forour parameters the model predicts a boundary that becomes flat at the leftand right, with a negative slope in between. This is borne out in Fig. 5.8(b)by the boundaries at Vm = 0.075 and 0.025 µm/s; in each case, the smalland large strain limits agree with the above estimations to within 1%.5.4.4 Comparison with experimentsThe critical condition for fluidization has been a focus of experimen-tal investigations, and some of the data can be compared with our modelpredictions. So far, disassembly of stress fibers has been investigated un-der several type of mechanical perturbations: compress-hold (CH) [23, 182],stretch-compress (SC)[18, 59, 109], compress-stretch (CS) [18], and cyclicstretch [89, 96–99, 212]. Given sufficiently severe deformation (in terms ofε0 or ε˙0), all these maneuvers have been shown to induce disassembly ofstress fibers in the direction of compression.Figure 5.8(a) compares the predicted critical condition for the SC ma-neuver with the experimental data of Chen et al. [18]. Since the cyclicstretch amounts to repeated cycles of the SC maneuver, its critical (ε0, ε˙0)map should be essentially the same as that for the SC maneuver. Thus, wehave included the cyclic-stretching data of Tondon et al. [212] and Kaunaset al. [99] as well. In these two data sets, we take the onset of SF alignmentin the orthogonal direction as indicative of stress fiber disassembly in thestretching direction. First, note that all the data are for relatively smallstrain, ε˙ ≤ 0.1, and thus cannot validate our prediction of a positive slopefor larger strains. Within this range, the data do suggest a boundary cor-responding to the critical condition that has a negative slope, as predictedby the model. Quantitatively, the predicted boundary is consistent with theonly data point for SC [18], but is considerably higher than that suggestedby the cyclic-stretch data.Now we turn to the CS and CH maneuvers in Fig. 5.8(b). Again, theexperimental data of Chen et al. [18] is consistent with the predicted criticalcondition; no fluidization occurs for ε0 = 0.1 and ε˙ = 0.05 s−1, a conditionbelow the theoretical curve. However, the CH experiments of Costa et al.[23] and Sato et al. [182] have reported buckling or disassembly of SFs atlower strain rates of 0.02–0.03 s−1. The data of Costa et al. [23] suggest aboundary with a negative slope, as predicted by the model, but the experi-mental boundary is much below the model prediction.925.5. ConclusionComparing the experiments cited in Fig. 5.8 and similar ones in the liter-ature, one may divide them into two groups according to substrate rigidity:experiments by Fredberg and colleagues [18, 59, 109, 213] have used verysoft substrates, typically collagen and polyacrylamide gels, with a Young’smodulus E up to 4 kPa, while all the other experiments [23, 99, 182, 212]have used much stiffer silicone-rubber substrates with E on the order ofMPa. It is well known that substrate rigidity affects the cytoskeletal struc-ture and mechanical behavior of adherent cells [37, 68]. As we have usedVm = 0.075 µm/s in the model, a value matched to the excess tension inthe soft-substrate experiment of Trepat et al. [213] (Fig. 5.7), it is reason-able that the predicted critical condition should agree with that of Chenet al. [18], another soft-substrate experiment. Stiffer substrates induce alarger number of stress fibers in the cell, larger traction forces and greatermodulus of the whole cell [92, 207]. In turn, these should affect the myosinconcentration and sliding speed. Unfortunately, it is difficult to representsuch effects in the model, say through Eq. (5.6), with any confidence. As-suming that the stiff substrates [23, 99, 182, 212] engender a stiffer myosinapparatus and slower myosin relaxation, we have computed the critical con-ditions for fluidization for a smaller characteristic Vm = 0.025 µm/s. Plottedas a dashed line in Fig. 5.8(b), the slower myosin relaxation does bring themodel prediction closer to these stiff-substrate experiments for the CS andCH maneuvers.To sum up the comparison above, the model prediction agrees qualita-tively, and in some cases semi-quantitatively, with experimentally measuredcritical strains and strain rates for SF disassembly. Encouraging as it is,the agreement should be considered preliminary and perhaps tentative, inview of the uncertainty in the parameter Vm as well as the various assump-tions that have gone into the model. These need to be investigated morethoroughly in future work.5.5 ConclusionThe response of cells to external strain is a complex process. Prior ex-periments have explored its various aspects, and the results are not always inagreement. Therefore, we have chosen a well-defined feature of the processfor theoretical modeling: the sensitivity of cell fluidization to strain history.More specifically, cells fluidize after a stretch-compress (SC) maneuver butnot after a compress-stretch (CS) maneuver of the same strain magnitudeand strain rate. Logically, our model is based on the following chain of935.5. Conclusionreasoning, each step being supported by experimental evidence explainedin the preceding sections. First, cell fluidization is due to disassembly ofthe cytoskeletal network, and in particular the stress fibers. Second, thestress fibers disintegrate as a result of the detachment of myosin motors enmasse from the actin filaments. Third, the detachment of myosin motors isdue to loss of tensile force in the stress fibers. Finally, the loss of tensionand appearance of a compressive force stem from the external strain beingimposed on the cell. To elucidate the dependence of fluidization on strainhistory, we have focused on modeling a single stress fiber undergoing thevarious protocols of stretching and compression. The model does not ac-count for the actual disintegration of the cytoskeleton, nor the subsequentresolidification. Chemical signaling has been neglected.Coupling the myosin crossbridge with a viscoelastic Kelvin-Voigt ele-ment, the model demonstrates a striking contrast between the SC and CSmaneuvers. In SC, the initial stretching leaves most of the myosin motorsin a post-power-stroke attached state, which rapidly disengage at the startof the compression stage, once the tension is lost. This is taken to be theprecursor of fluidization. In CS of the same magnitude, on the other hand,the initial compression has a much milder effect on the myosin population,and no fluidization results. We further explore the critical strain and strainrate for fluidization in SC, CS and compress-hold (CH) maneuvers and com-pare the model prediction with experiments. There is qualitative and evensemi-quantitative agreement. The main contribution of this work, therefore,is a mechanical explanation of how adherent cells react differently accordingto the history of the external strain.The Kelvin-Voigt-Myosin model proposed here is arguably the simplestmodel that allows one to describe the viscosity, elasticity and contractilityof a stress fiber, and to capture the mechanism of strain-history dependence.For this reason, it must be viewed as only one ingredient in a more compre-hensive description of the process of fluidization and resolidification. Severalimportant factors are left out. For example, we take the loss of tension asmarking the onset of fluidization. The subsequent events of stress-fiber buck-ling and disassembly and depolymerization of the actin filaments are notexplicitly accounted for. Moreover, resolidification—the longer-term recov-ery of cell stiffness—presents interesting mechanical and biological questionsthat we have not attempted to address. Experiments suggest that the slowrecovery of filamentous actin and actin-myosin connectivity is ATP-drivenand regulated by signal pathways, with zyxin being a potentially key player[18, 109]. Due to the model’s simplicity, it cannot capture the time scalefor the relaxation of cell modulus after a stretch-hold maneuver, nor can it945.5. Conclusionreflect the effect of substrate stiffness, which influences the mechanosensingand cytoskeletal remodeling of adherent cells. Finally, the model ignoresa potentially important mechanism for redistributing myosin motors fromthe cytoplasm to the cytoskeleton under external forcing [122, 123, 198]. Inthis context, the current model should be seen as a preliminary effort atformulating and rationalizing the response of cells to mechanical forcing.95Chapter 6Conclusions andrecommendationsIn this chapter, I summary the major findings and contributions of myresearch. I also list limitations of my current work and give recommendationsfor future work.6.1 Malaria-infected red blood cellsIn this project, I coupled a discrete mechanical cell model with a paral-lel SPH solver. The model was applied to simulate the traverse of infectedred blood cells (iRBCs) through microfluidic channels. In this study, I in-vestigated how stiffening of the membrane, presence of parasite cluster andmorphological changes of the iRBCs affect the cell passing through the mi-crofluidic channels. The numerical results revealed that the presence of therigid parasite cluster plays the major role in causing the blockage of mi-crofluidic flows. The morphological changes of the iRBCs also significantlyreduce the iRBCs’ deformability.This study suffers from some limitations. For example, the model doesnot account for adhesive properties of the iRBC membrane. The viscosity ofthe membrane has been neglected. The cytoplasmic viscosity was treated asthe same as that of the surrounding fluid. The spectrin-based cytoskeletonunderlying the RBC membrane is treated as an elastic triangular network,and its dynamics and potential remodeling have not been accounted in themodel. Regarding the methodology, the SPH method is computationallycostly due to its explicit time advancing scheme. Especially for the problemwith small length scale, the major time constraint comes from the viscousdiffusion.Regarding the limitations reviewed above, there is still much room for usto improve the current cell model. The most interesting improvement willbe to include the dynamics of the spectrin network. In reality, the red cellcytoskeleton undergoes structural remodeling [116, 243], like binding and966.2. Activation of neutrophilsunbinding between the spectrins. This property contributes to the RBCdeformability. Thanks to the discrete nature of the model, it will be rel-atively straightforward to include such a restructuring mechanism for thecytoskeleton.6.2 Activation of neutrophilsIn this project, I investigated the activation of the neutrophils whilepassing through a microfluidic channel. The numerical results show thatthe activation of the neutrophils relies on two factors: (1) mechanical defor-mation disassembling the cytoskeleton, which removes the cortical tensionthat would inhibit local actin buildup, and (2) sufficiently slow entry ofthe cell into the channel, which allows the chemical kinetics to develop thesignal polarity. This provides what appears to be the first comprehensiveexplanation of a set of experimental observations.The model still suffers from some limitations. Only one signal proteinRac has been used to identify the signal polarization. The chemical kineticswere only accounted for on the cell membrane. The model simply treatedthe cortical tension and the protrusion force as functions of the Rac levelon the membrane, thus overlooking the detailed pathways leading from Racto actomyosin buildup. Limited by the 2D simulation, the rate of the mem-brane bending was used as the disassembling criterion for the cytoskeleton.Actual mechanics of the process is disregarded.The cell may respond distinctively to different mechanical deformations.This is because certain mechanosensing proteins can be activated by a spe-cific type of deformation [122]. For example, filamin A (FLNa) is sensitiveto shear deformation. In the future, the 2D model can be extended to a 3Done. By doing so, the cell response to different deformations can be inves-tigated. In the current model, the cortical tension inhibition is implicitlyaccounted for by a coarse-grained parameter kτ . In the future, it would beinteresting to simulate the FLNa-FilGAP pathway directly.6.3 Cell fluidization and reinforcementIn this project, I coupled the cross-bridge cycle of myosin motors witha viscoelastic Kelvin-Voigt element to represent the stress fiber. The dy-namics of the model under a stretch-compress (SC) and compress-stretch(CS) deformation with 10% magnitude were investigated. The model suc-cessfully predicted that under the SC maneuver the stretch increases the976.3. Cell fluidization and reinforcementstress-fiber (SF) tension and suppresses the myosin detachment. The subse-quent compression then causes a large proportion of the myosin populationto disengage rapidly from actin filaments.Though the model captured the asymmetric behavior of the SF underthe SC and CS maneuvers, it is still a drastically simplified representation ofthe cell-stretching experiments. As such, it does not explicitly account forthe disassembly and assembly of the SFs, but uses the loss of tension as anindication of the onset of disassembly. Besides, it does not account for thebiochemical signal interactions of the focal adhesions and the stress fibers.For future work, the individual SF model can be extended to multi-SF models. The reorientation of the SFs under cyclic stretches has beenobserved experimentally [99, 108]. The reorientation of the SFs is crucial tothe reorientation of the adherent cells under cyclic stretching. It involves thedynamics of the focal adhesions and the disassembly/reassembly of the SFs.The single KVM model used here displays complex responses to mechanicaldeformations. It would be interesting to couple the KVM model with adiscrete particle-spring system [12, 219, 236] to represent bundles of SFsand even the cytoskeleton. Based on that, the SF reorientation dynamicscan be investigated on a much more detailed level.98Bibliography[1] Manouk Abkarian and Annie Viallat. Vesicles and red blood cells inshear flow. Soft Matter, 4:653–657, 2008.[2] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter.DNA, Chromosomes, and Genomes. In Molecular Biology of the Cell,chapter 4, pages 191–234. Garland Scientific, 2002.[3] J. Alvarado, M. Sheinman, A. Sharma, F. C. MacKintosh, and G. H.Koenderink. Molecular motors robustly drive active gels to a criticallyconnected state. Nat. Phys., 9:591–597, 2013.[4] M. Antia, T. Herricks, and P. K. Rathod. Microfluidic approaches tomalaria pathogenesis. Cell. Microbiol., 10:1968–1974, 2008.[5] F. P. T. Baaijens, W. R. Trickey, T. A. Laursen, and F. Guilak. Largedeformation finite element analysis of micropipette aspiration to de-termine the mechanical properties of the chondrocyte. Ann. Biomed.Eng., 33:491–501, 2005.[6] D. Barthe`s-Biesel, A. Diaz, and E. Dhenin. Effect of constitutive lawsfor two-dimensional membranes on flow-induced capsule deformation.J. Fluid. Mech., 460:211–222, 2002.[7] J. Bates, S. Bullimore, A. Politi, J. Sneyd, R. Anafi, and A. Lauzon.Transient oscillatory force-length behavior of activated airway smoothmuscle. Am. J. Physiol.-Lung C., 297:L362–L372, 2009.[8] A. Besser and U. S. Schwarz. Coupling biochemistry and mechanicsin cell adhesion: a model for inhomogeneous stress fiber contraction.New J. Phys., 9:425, 2007.[9] D. H. Boal and M. Rao. Topology changes in fluid membranes. Phys.Rev. A, 46:3037–3045, 1992.99Bibliography[10] G. Boedec, M. Leonetti, and M. Jaeger. 3D vesicle dynamics simula-tions with a linearly triangulated surface. J. Comput. Phys., 230:1020–1034, 2011.[11] S. K. Boey, D. H. Boal, and D. E. Discher. Simulations of the erythro-cyte cytoskeleton at large deformation. I. Microscopic models. Bio-phys. J., 75:1573–1583, 1998.[12] C. Borau, T. Kim, T. Bidone, J. M. Garc´ıa-Aznar, and R. D. Kamm.Dynamic mechanisms of cell rigidity sensing: Insights from a compu-tational model of actomyosin networks. PLoS One, 7:e49174, 2012.[13] H. Bow, I. V. Pivkin, M. Diez Silva, S. J. Goldfless, M. Dao, J. C.Niles, S. Suresh, and J. Han. A microfabricated deformability-basedflow cytometer with application to malaria. Lab Chip, 11:1065–1073,2011.[14] B. S. Brook. Emergence of airway smooth muscle mechanical behaviorthrough dynamic reorganization of contractile units and force trans-mission pathways. J. Appl. Physiol., 116:980–997, 2014.[15] R. A. Brown, R. Prajapati, D. A. McGrouther, I. V. Yannas, andM. Eastwood. Tensional homeostasis in dermal fibroblasts: mechanicalresponses to mechanical loading in three-dimensional substrates. J.Cell. Physiol., 175:323–332, 1998.[16] P. A. Buffet, I. Safeukui, G. Deplaine, V. Brousse, V. Prendki,M. Thellier, G. D. Turner, and O. M. Puijalon. The pathogenesisof Plasmodium falciparum malaria in humans: insights from splenicphysiology. Blood, 117:381–392, 2011.[17] B. Chen, R. Kemkemer, M. Deibler, J. Spatz, and H. Gao. Cyclicstretch induces cell reorientation on substrates by destabilizing catchbonds in focal adhesions. PLoS One, 7:e48346, 2012.[18] C. Chen, R. Krishnan, E. Zhou, A. Ramachandran, D. Tambe, K. Ra-jendran, R. M. Adam, L. Deng, and J. J. Fredberg. Fluidization andresolidification of the human bladder smooth muscle cell in responseto transient stretch. PLoS One, 5:e12035, 2010.[19] L. Chin, P. Yue, J. J. Feng, and C. Y. Seow. Mathematical simulationof muscle crossbridge cycle and force-velocity relationship. Biophys.J., 91:3653–3663, 2006.100Bibliography[20] G. R. Cokelet and H. J. Meiselman. Rheological comparison ofhemoglobin solutions and erythrocyte suspensions. Science, 162:275–277, 1968.[21] J. Colombelli, A. Besser, H. Kress, E. G. Reynaud, P. Girard,E. Caussinus, U. Haselmann, J. V. Small, U. S. Schwarz, and E. H.Stelzer. Mechanosensing in actin stress fibers revealed by a close cor-relation between force and protein localization. J. Cell Sci., 122:1665–1679, 2009.[22] M. A. Conti, S. Kawamoto, and R. S. Adelstein. Non-muscle myosinII. 2008. In Proteins and Cell Regulation, Vol. 7. A. J. Ridley and J.Frampton, editors. Springer, Berlin. 223–264.[23] K. D. Costa, W. J. Hucker, and F. C. P. Yin. Buckling of actin stressfibers: a new wrinkle in the cytoskeletal tapestry. Cell Motil. Cytoskel.,52:266–274, 2002.[24] M. F. Coughlin and G. W. Schmid-Scho¨nbein. Pseudopod projec-tion and cell spreading of passive leukocytes in response to fluid shearstress. Biophys. J., 87:2035–2042, 2004.[25] G. Coupier, B. Kaoui, T. Podgorski, and C. Misbah. Noninertiallateral migration of vesicles in bounded poiseuille flow. Phys. Fluids,20:111702, 2008.[26] A. E. Cowan, L. Nakhimovsky, D. G. Myles, and D. E. Koppel. Barri-ers to diffusion of plasma membrane proteins form early during guineapig spermiogenesi. Biophys. J., 73:507–516, 1997.[27] H. A. Cranston, C. W. Boylan, G. L. Carroll, S. P. Sutera, J. R.Williamson, I. Y. Gluzman, and D. J. Krogstad. Plasmodium falci-parum maturation abolishes physiologic red cell deformability. Sci-ence, 223:400–403, 1984.[28] M. Dao, J. Li, and S. Suresh. Molecularly based analysis of deforma-tion of spectrin network and human erythrocyte. Mater. Sci. Eng. C,26:1232–1244, 2006.[29] M. Dao, C. T. Lim, and S. Suresh. Mechanics of the human red bloodcell deformed by optical tweezers. J. Mech. Phys. Solids, 51:2259 –2280, 2003.101Bibliography[30] R. De, A. Zemel, and S. A. Safran. Dynamics of cell orientation. Nat.Phys., 3:655–659, 2007.[31] E. P. Debold, J. B. Patlak, and D. M. Warshaw. Slip sliding away:Load-dependence of velocity generated by skeletal muscle myosinmolecules in the laser trap. Biophys. J., 89:L34–L36, 2005.[32] S. Deguchi, T. Ohashi, and M. Sato. Tensile properties of singlestress fibers isolated from cultured vascular smooth muscle cells. J.Biomech., 39:2603–2610, 2006.[33] V. S. Deshpande, R. M. McMeeking, and A. G. Evans. A bio-chemo-mechanical model for cell contractility. PNAS, 103:14015–14020, 2006.[34] V. S. Deshpande, R. M. McMeeking, and A. G. Evans. A model forthe contractility of the cytoskeleton including the effects of stress-fibreformation and dissociation. Proc. R. Soc., 463:787–815, 2007.[35] V. S. Deshpande, M. Mrksich, R. M. McMeeking, and A. G. Evans. Abio-mechanical model for coupling cell contractility with focal adhesionformation. J. Mech. Phys. Solids, 56:1484–1510, 2008.[36] D. E. Discher, D. H. Boal, and S. K. Boey. Simulation of the erythro-cyte cytoskeleton at large deformation. II. Micropipette aspiration.Biophys. J., 75:1584–1597, 1998.[37] D. E. Discher, P. Janmey, and Y. Wang. Tissue cells feel and respondto the stiffness of their substrate. Science, 310:1139–1143, 2005.[38] A. M. Dondorp, E. Pongponratn, and N. J. White. Reduced micro-circulatory flow in severe falciparum malaria: pathophysiology andelectron-microscopic pathology. Acta Trop., 89:309–317, 2004.[39] C. Dong and R. Skalak. Leukocyte deformability: Finite element mod-eling of large viscoelastic deformation. J. Theor. Biol., 158:173–193,1992.[40] C. Dong, R. Skalak, K. L. Sung, G. W. Schmid-Scho¨nbein, and S. C.Chien. Passive deformation analysis of human leukocytes. J. Biomech.Eng., 110:27–36, 1988.[41] G. M. Donovan. Modelling airway smooth muscle passive length adap-tation via thick filament length distributions. J. Theor. Biol., 333:102–108, 2013.102Bibliography[42] J. L. Drury and M. Dembo. Hydrodynamics of micropipette aspira-tion. Biophys. J., 78:110–128, 1999.[43] J. L. Drury and M. Dembo. Aspiration of human neutrophils: Effectsof shear thinning and cortical dissipation. Biophys. J., 81:3166–3177,2001.[44] Jules Dupire, Marius Socol, and Annie Viallat. Full dynamics of a redblood cell in shear flow. PNAS, 109:20808–20813, 2012.[45] T. T. Egelhoff, T. V. Naismith, and F. V. Brozovich. Myosin-basedcortical tension in dictyostelium resolved into heavy and light chain-regulated components. J. Muscle Res. Cell Motil., 17:269–274, 1996.[46] A. J. Ehrlicher, F. Nakamura, J. H. Hartwig, D. A. Weitz, and T. P.Stossel. Mechanical strain in actin networks regulates FilGAP andintegrin binding to filamin A. Nature, 478:260–263, 2011.[47] A. Esposito, J. Choimet, J. Skepper, J. Mauritz, V. Lew, C. Kaminski,and T. Tiffert. Quantitative imaging of human red blood cells infectedwith Plasmodium falciparum. Biophys. J., 99:953–960, 2010.[48] A. Esposito, T. Tiffert, J. M. A. Mauritz, S. Schlachter, L. H. Ban-nister, C. F. Kaminski, and V. L. Lew. FRET imaging of hemoglobinconcentration in Plasmodium falciparum-infected red cells. PLoS One,3:e3780, 2008.[49] E. Evans. New membrane concept applied to the analysis of fluid shear- and micropipette-deformed red blood cells. Biophys. J., 13:941–954,1973.[50] E. Evans and Y. Fung. Improved measurements of the erythrocytegeometry. Microvasc. Res., 4:335–347, 1972.[51] E. Evans and A. Yeung. Apparent viscosity and cortical tension ofblood granulocytes determined by micropipette aspiration. Biophys.J., 56:151–160, 1989.[52] E. A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. CombinedImmersed-Boundary finite-difference methods for three-dimensionalcomplex flow simulations. J. Comput. Phys., 161:35–60, 2000.[53] D. A. Fedosov, B. Caswell, and G. E. Karniadakis. A multiscale redblood cell model with accurate mechanics, rheology, and dynamics.Biophys. J., 98:2215–2225, 2010.103Bibliography[54] D. A. Fedosov, B. Caswell, and G. E. Karniadakis. Wall shear stress-based model for adhesive dynamics of red blood cells in malaria. Bio-phys. J., 100:2084–2093, 2011.[55] D. A. Fedosov, B. Caswell, S. Suresh, and G. E. Karniadakis.Quantifying the biophysical characteristics of Plasmodium-falciparum-parasitized red blood cells in microcirculation. PNAS, 108:35–39, 2011.[56] T. M. Fischer, M. Stohr-Liesen, and H. Schmid-Schonbein. The red cellas a fluid droplet: Tank tread-like motion of the human erythrocytemembrane in shear flow. Science, 202:894–896, 1978.[57] D. A. Fletcher and R. D. Mullins. Cell mechanics and the cytoskeleton.Nature, 463:485–492, 2010.[58] S. Gabriele, A. M. Benoliel, P. Bongrand, and O. The´odoly. Mi-crofluidic investigation reveals distinct roles for actin cytoskeletonand myosin ii activity in capillary leukocyte trafficking. Biophys. J.,96:4308–4318, 2009.[59] N. Gavara, P. Roca-Cusachs, R. Sunyer, R. Farre´, and D. Navajas.Mapping cell-matrix stresses during stretch reveals inelastic reorgani-zation of the cytoskeleton. Biophys. J., 95:464–471, 2008.[60] G. Giannone, G. Jiang, D. H. Sutton, D. R. Critchley, and M. P.Sheetz. Talin1 is critical for force-dependent reinforcement of initialintegrin cytoskeleton bonds but not tyrosine kinase activation. J. CellBiol., 163:409–419, 2003.[61] S. C. Gifford, M. G. Frank, J. Derganc, C. Gabel, R. H. Austin,T. Yoshida, and M. W. Bitensky. Parallel microchannel-based mea-surements of individual erythrocyte areas and volumes. Biophys. J.,84:623–633, 2003.[62] F. K. Glenister, R. L. Coppel, A. F. Cowman, N. Mohandas, andB. M. Cooke. Contribution of parasite proteins to altered mechanicalproperties of malaria-infected red blood cells. Blood, 99:1060–1063,2002.[63] F. K. Glenister, K. M. Fernandez, L. M. Kats, E. Hanssen, N. Mo-handas, R. L. Coppel, and B. M. Cooke. Functional alteration of redblood cells by a megadalton protein of Plasmodium falciparum. Blood,113:919–928, 2009.104Bibliography[64] M. Go´mez-Gesteira, A. J. C. Crespo, B. D. Rogers, R. A. Dalrym-ple, J. M. Dominguez, and A. Barreiro. Sphysics – development of afree-surface fluid solver – Part 2: Efficiency and test cases. Comput.Geosci., 48:300–307, 2012.[65] M. Go´mez-Gesteira, B. D. Rogers, A. J. C. Crespo, R. A. Dalrymple,M. Narayanaswamy, and J. M. Dominguez. Sphysics – development ofa free-surface fluid solver – Part 1: Theory and formulations. Comput.Geosci., 48:289–299, 2012.[66] Q. Guo, S. P. Duffy, K. Matthews, A. T. Santoso, M. D. Scott,and H. Ma. Microfluidic analysis of red blood cell deformability. J.Biomech., 47:1767–1776, 2014.[67] Q. Guo, S. J. Reiling, P. Rohrbach, and H. Ma. Microfluidic biome-chanical assay for red blood cells parasitized by Plasmodium falci-parum. Lab Chip, 12:1143–1150, 2012.[68] W.-H. Guo, Margo T. Frey, Nancy A. Burnham, and Y.-L. Wang.Substrate rigidity regulates the formation and maintenance of tissues.Biophys. J., 90:2213–2220, 2006.[69] C. M. Hai and R. A. Murphy. Cross-bridge phosphorylation and regu-lation of latch state in smooth muscle. Am. J. Physiol., 254:C99–C106,1988.[70] A. Hall. Rho GTPases and the actin cytoskeleton. Science, 279:509–513, 1998.[71] E. Hanssen, P. J. McMillan, and L. Tilley. Cellular architectureof Plasmodium falciparum-infected erythrocytes. Int. J. Parasitol.,40:1127–1135, 2010.[72] W. Helfrich. Elastic properties of lipid bilayers: theory and possibleexperiments. Z. Naturforsch. C, 28:693–703, 1973.[73] M. Herant, V. Heinrich, and M. Dembo. Mechanics of neutrophilphagocytosis: behavior of the cortical tension. J. Cell Sci., 118:1789–1797, 2005.[74] M. Herant, V. Heinrich, and M. Dembo. Mechanics of neutrophilphagocytosis: experiments and quantitative models. J. Cell Sci.,119:1903–1913, 2006.105Bibliography[75] M. Herant, W. A. Marganski, and M. Dembo. The mechanics ofneutrophils: Synthetic modeling of three experiments. Biophys. J.,84:3389–3413, 2003.[76] T. Herricks, M. Antia, and P. K. Rathod. Deformability limits of Plas-modium falciparum-infected red blood cells. Cell. Microbiol., 11:1340–1353, 2009.[77] T. Herricks, K. B. Seydel, M. Molyneux, T. Tayler, and P. K.Rathod. Estimating physical splenic filtration of Plasmodium falci-parum-infected red blood cell in malaria patients. Cell. Microbiol.,14:1880–1891, 2012.[78] A. V. Hill. The heat of shortening and the dynamic constants ofmuscle. Proc. Roy. Soc. Lond., 126:136–195, 1938.[79] H. Hinghofer-Szalkay and J. E. Greenlea. Continuous monitoring ofblood volume changes in humans. J. Appl. Physiol., 63:1003–1007,1987.[80] R. M. Hochmuth. Micropipette aspiration of living cells. J. Biomech.,33:15–22, 2000.[81] R. M. Hochmuth and R. E. Waugh. Erythrocyte membrane elasticityand viscosity. Ann. Rev. Physiol., 49:209–219, 1987.[82] B. D. Hoffman, C. Grashoff, and M. A. Schwartz. Dynamic molecularprocesses mediate cellular mechanotransduction. Nature, 475:316–323,2011.[83] W. R. Holmes, A. E. Carlsson, and L. Edelstein-Keshet. Regimesof wave type patterning driven by refractory actin feedback: transi-tion from static polarization to dynamic wave behavior. Phys. Biol.,9:046005, 2012.[84] S. M. Hosseini. A particle-based model for computing fluid flows andcell dynamics. PhD thesis, University of British Columbia, 2012.[85] S. M. Hosseini and J. J. Feng. A particle-based model for the transportof erythrocytes in capillaries. Chem. Eng. Sci., 64:4488–4497, 2009.[86] S. M. Hosseini and J. J. Feng. Pressure boundary conditions for com-puting incompressible flows with SPH. J. Comput. Phys., 230:7473–7487, 2011.106Bibliography[87] S. M. Hosseini and J. J. Feng. How malaria parasites reduce thedeformability of infected red blood cells. Biophys. J., 103:1–10, 2012.[88] H. W. Hou, W. C. Lee, M. C. Leong, S. Sonam, S. R. K. Vedula,and C. T. Lim. Microfluidics for applications in cell mechanics andmechanobiology. Cell. Mol. Bioeng., 4:591–602, 2011.[89] H. J. Hsu, C. F. Lee, A. Locke, S. Q. Vanderzyl, and R. Kaunas.Stretch-induced stress fiber remodeling and the activations JNK andERK depend on mechanical strain rate, but not FAK. PLoS One,5:e12470, 2010.[90] W. R. Dodson III and P. Dimitrakopoulos. Tank-treading of erythro-cytes in strong shear flows via a nonstiff cytoskeleton-based continuumcomputational modeling. Biophys. J., 99:2906–2916, 2010.[91] Y. Imai, H. Kondo, T. Ishikawa, C. T. Lim, and T. Yamaguchi. Mod-eling of hemodynamics arising from malaria infection. J. Biomech.,43:1386–1393, 2010.[92] R. A. Jannat, M. Dembo, and D. A. Hammer. Traction forces ofneutrophils migrating on compliant substrates. Biophys. J., 101:575–584, 2011.[93] G. R. Johnson and S. R. Beissel. Normalized smoothing function forSPH impact computations. Int. J. Numer. Methods Eng., 39:2725–2741, 1996.[94] K. Katoh, Y. Kano, M. Amano, H. Onishi, K. Kaibuchi, and K. Fu-jiwara. Rho-kinase mediated contraction of isolated stress fibers. J.Cell Biol., 153:569–583, 2001.[95] K. Katoh, Y. Kano, M. Masuda, H. Onishi, and K. Fujiwara. Isolationand contraction of the stress fiber. Mol. Biol. Cell, 9:1919–1938, 1998.[96] R. Kaunas and S. Deguchi. Multiple roles for myosin II in tensionalhomeostasis under mechanical loading. Cell. Mol. Bioeng., 4:182–191,2011.[97] R. Kaunas and H.-J. Hsu. A kinematic model of stretch-induced stressfiber turnover and reorientation. J. Theor. Biol., 257:320–330, 2009.[98] R. Kaunas, H.-J. Hsu, and S. Deguchi. Sarcomeric model of stretch-induced stress fiber reorganization. Cell Health Cytoskel., 3:13–22,2011.107Bibliography[99] R. Kaunas, P. Nguyen, S. Usami, and S. Chien. Cooperative effectsof rho and mechanical stretch on stress fiber organization. Proc. Natl.Acad. Sci. USA, 102:15895–15900, 2005.[100] G. Ke´sma´rky, P. Kenyeres, M. Ra´bai, and K. To´th. Plasma viscosity:a forgotten variable. Clin. Hemorheol. Micro., 39:243–246, 2008.[101] T. Kim, W. Hwang, and R. D. Kamm. Computational analysis of across-linked actin-like networks. Exp. Mech., 49:91–104, 2009.[102] T. Kim, W. Hwang, H. Lee, and R. D. Kamm. Computational analysisof viscoelastic properties of crosslinked actin networks. PLoS Comput.Biol., 5:e1000439, 2009.[103] K. Kimura, M. Ito, M. Amano, K. Chihara, Y. Fukata, M. Naka-fuku, B. Yamamori, J. Feng, T. Nakano, K. Okawa, A. Iwamatsu, andK. Kaibuchi. Regulation of myosin phosphatase by Rho and Rho-associated kinase (Rho-kinase). Science, 273:245–248, 1996.[104] E. Kolaczkowska and P. Kubes. Neutrophil recruitment and functionin health and inflammation. Nat. Rev. Immunol., 13:159–175, 2013.[105] J. Kolega. Cytoplasmic dynamics of myosin ii a and ii b: spatialsorting of isoforms in locomoting cells. J. Cell Sci., 111:2085–2095,1998.[106] J. Kolega. Asymmetric distribution of myosin iib in migration en-dothelial cells is regulated by a rho-dependent kinase and contributesto tail retraction,. Mol. Biol. Cell, 14:4745–4757, 2003.[107] M. Kova´cs, K. Thirumurugan, P. J. Knight, and J. R. Sellers. Load-dependent mechanism of nonmuscle myosin 2. Proc. Natl. Acad. Sci.USA, 104:9994–9999, 2007.[108] R. Krishnan, E. P. Canovic´, A. L. Iordan, K. Rajendran, G. Manomo-han, A. P. Pirentis, M. L. Smith, J. P. Butler, J. J. Fredberg, andD. Stamenovic´. Fluidization, resolidification, and reorientation of theendothelial cell in response to slow tidal stretches. Am. J. Physiol.Cell Physiol., 303(4):C368–C375, 2012.[109] R. Krishnan, C. Y. Park, Y. C. Lin, J. Mead, R. T. Jaspers, X. Trepat,G. Lenormand, D. Tambe, A. V. Smolensky, A. H. Knoll, J. P. Butler,and J. J. Fredberg. Reinforcement versus fluidization in cytoskeletalmechanoresponsiveness. PLoS One, 4:e5486, 2009.108Bibliography[110] S. Kumar, I. Z. Maxwell, A. Heisterkamp, T. R. Polte, T. P. Lele,M. Salanga, E. Mazur, and D. E. Ingber. Viscoelastic retraction ofsingle living stress fibers and its impact on cell shape, cytoskeletal or-ganization, and extracellular matrix mechanics. Biophys. J., 90:3762–3773, 2006.[111] E. Lac, D. Barthe`s-Biesel, N. A. Pelekasis, and J. Tsamopoulos. Spher-ical capsules in threedimensional unbounded stokes flow: effect of themembrane constitutive law and onset of buckling. J. Fluid. Mech.,516:303–334, 2004.[112] E. Lac, A. Morel, and D. Barthe`s-Biesel. Hydrodynamic interactionbetween two identical capsules in simple shear flow. J. Fluid. Mech.,573:149–169, 2007.[113] E. S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, andP. Stansby. Comparisons of weakly compressible and truly incompress-ible algorithms for the SPH mesh free particle method. J. Comput.Phys., 227:8417–8436, 2008.[114] G. Y. H. Lee and C. T. Lim. Bimechanics approaches to studyinghuman diseases. Trends Biotechnol., 25:111–118, 2007.[115] F. Y. Leong, Q. Li, C. T. Lim, and K. H. Chiam. Modeling cell entryinto a micro-channel. Biomech. Model Mechanobiol, 10:755–766, 2011.[116] J. Li, G. Lykotrafitis, M. Dao, and S. Suresh. Cytoskeletal dynamicsof human erythrocyte. PNAS, 104:4937–4942, 2007.[117] L. D. Libersky, A. G. Petschek, T. C. Carney, J. R. Hipp, and F. A.Allahdadi. High strain lagrangian hydrodynamics: a three dimensionalsph code for dynamic material response. J. Comput. Phys., 109:67–75,1993.[118] C. T. Lim, E. H. Zhou, and S. T. Quek. Mechanical models for livingcells—a review. J. Biomech., 39:195–216, 2006.[119] B. Lin, W. R. Holmes, C. J. Wang, T. Ueno, A. Harwell, L. Edelstein-Keshet, T. Inoue, and A. Levchenko. Synthetic spatially graded Racactivation drives cell polarization and movement. PNAS, 26:E3668–E3677, 2012.[120] G. R. Liu and M. B. Liu. Smoothed particle hydrodynamics: a meshfreeparticle method. World Scientific Publishing Co. Pte. Ltd., 2003.109Bibliography[121] L. Lu, Y. Feng, W. J. Hucker, S. J. Oswald, G. D. Longmore, andF. C. P. Yin. Actin stress fiber pre-extension in human aortic en-dothelial cells. Cell Motil. Cytoskel., 65:281–294, 2008.[122] T. Luo, K. Mohan, P. A. Iglesias, and D. N. Robinson. Molecularmechanisms of cellular mechanosensing. Nat. Mater., 12:1064–1071,2013.[123] T. Luo, K. Mohan, V. Srivastava, Y. Ren, P. A. Iglesias, andD. N. Robinson. Understanding the cooperative interaction betweenmyosin ii and actin cross-linkers mediated by actin filaments duringmechanosensation. Biophys. J., 102:238–247, 2012.[124] L. M. Machesky and A. Hall. Role of actin polymerization and ad-hesion to extracellular matrix in Rac- and Rho-induced cytoskeletalreorganization. J. Cell Biol., 138:913–926, 1997.[125] A. G. Maier, B. M. Cooke, A. F. Cowman, and L. Tilley. Malaria par-asite proteins that remodel the host erythrocyte. Nat. Rev. Microbiol.,7:341–354, 2009.[126] A. Makino, M. Glogauer, G. M. Bokoch, S. Chien, and G. W. Schmid-Scho¨nbein. Control of neutrophil pseudopods by fluid shear: role ofRho family GTPases. Am. J. Physiol. Cell Physiol., 288:C863–C871,2005.[127] A. F. M. Mare´e, A. Jilkine, A. Dawes, V. A. Grieneisen, andL. Edelstein-Keshet. Polarization and movement of keratocytes: Amultiscale modelling approach. B. Math. Biol., 68:1169–1211, 2006.[128] S. V. Marella and H. S. Udaykumar. Computational analysis of thedeformability of leukocytes modeled with viscous and elastic structuralcomponents. Phys. Fluids, 16:244–264, 2004.[129] Bryan T. Marshall, Mian Long, James W. Piper, Tadayuki Yago,Rodger P. McEver, and Cheng Zhu. Direct observation of catch bondsinvolving cell-adhesion molecules. Nature, 423:190–194, 2003.[130] M. A. Mata, M. Dutot, L. Edelstein-Keshet, and W. R. Holmes. Amodel for intracellular actin waves explored by nonlinear local pertur-bation analysis. J. Theor. Biol., 334:149–161, 2013.110Bibliography[131] T. S. Matsui, K. Ito, R. Kaunas, M. Kanzaki, M. Sato, and S. Deguchi.Non-muscle myosin II induces disassembly of actin stress fibres inde-pendently of myosin light chain dephosphorylation. Interface Focus,1:754–766, 2011.[132] T. S. Matsui, K. Ito, R. Kaunas, M. Sato, and S. Deguchi. Actinstress fibers are at a tipping point between conventional shorteningand rapid disassembly at physiological levels of Mg-ATP. Biochem.Biophys. Res. Commun., 395:301–306, 2010.[133] B. D. Matthews, D. R. Overby, R. Mannix, and D. E. Ingber. Cellularadaptation to mechanical stress: role of integrins, Rho, cytoskeletaltension and mechanosensitive ion channels. J. Cell Sci., 119(3):508–518, 2006.[134] S. M. McFaul, B. K. Lin, and H. Ma. Cell separation based on size anddeformability using microfluidic funnel ratchets. Lab Chip, 12:2369–2376, 2012.[135] S. Mijailovich, J. Butler, and J. J. Fredberg. Perturbed equilibria ofmyosin binding in airway smooth muscle: bond-length distributions,mechanics, and ATP metabolism. Biophys. J., 79:2667–2681, 2000.[136] S. M. Mijailovich, M. Kojic, M. Zivkovic, B. Fabry, and J. J. Fred-berg. A finite element model of cell deformation during magnetic beadtwisting. J. Appl. Physiol., 93:1429–1436, 2002.[137] L. H. Miller, D. I. Baruch, K. Marsh, and O. K. Doumbo. Thepathogenic basis of malaria. Nature, 415:673–679, 2002.[138] R. Mittal and G. Laccarino. Immersed boundary method. Annu. Rev.Fluid Mech., 37:239–261, 2005.[139] T. Mizutani, H. Haga, and K. Kawabata. Cellular stiffness responseto external deformation: tensional homeostasis in a single fibroblast.Cell Motil. Cytoskel., 59:242–248, 2004.[140] F. Moazzam, F. A. Delano, B. W. Zweifach, and G. W. Schmid-Scho¨nbein. The leukocyte response to fluid stress. PNAS, 94:5338–5343, 1997.[141] N. Mohandas and E. Evans. Mechanical properties of the red cell mem-brane in relation to molecular structure and genetic defects. Annu.Rev. Biophys. Biomol. Struct., 23:787–818, 1994.111Bibliography[142] J. J. Monaghan. Smoothed particle hydrodynamics. Ann. Rev. Astron.Astrophys., 30:543–574, 1992.[143] J. J. Monaghan. Simulating free surface flows with SPH. J. Comput.Phys., 110:399–406, 1994.[144] J. J. Monaghan. Smoothed particle hydrodynamics. Rep. Prog. Phys.,68:1703–1759, 2005.[145] J. J. Monaghan. Smoothed particle hydrodynamics and its diverseapplications. Ann. Rev. Fluid Mech., 44:323–346, 2012.[146] Y. Mori, A. Jilkine, and L. Edelstein-Keshet. Wave-pinning andcell polarity from a bistable reaction-diffusion system. Biophys. J.,94:3684–3697, 2008.[147] Y. Mori, A. Jilkine, and L. Edelstein-Keshet. Asymptotic and bifur-cation analysis of wave-pinning in a reaction-diffusion model for cellpolarization. SIAM J. Appl. Math., 71:1401–1427, 2011.[148] N. Morone, T. Fujiwara, K. Murase, R. S. Kasai, H. Ike, S. Yuasa,J. Usukura, and A. Kusumi. Three-dimensional reconstruction of themembrane skeleton at the plasma membrane interface by electron to-mography. J. Cell Biol., 174:851–862, 2006.[149] K. I. Morozov and L. M. Pismen. Cytoskeleton fluidization versusresolidification: Prestress effect. Phys. Rev. E, 83:051920, 2011.[150] K. I. Morozov and L. M. Pismen. Strain dependence of cytoskeletonelasticity. Soft Matter, 8:9193–9199, 2012.[151] M. P. Murrell and M. L. Gardel. F-actin buckling coordinates con-tractility and severing in a biomimetic actomyosin cortex. PNAS,109:20820–20825, 2012.[152] F. Nakamura. FilGAP and its close relatives: a mediator of Rho-Racantagonism that regulates cell morphology and migration. Biochem.J., 453:17–25, 2013.[153] G. B. Nash, E. O’Brien, E. C. Gordon-Smith, and J. A. Dormandy.Abnormalities in the mechanical properties of red blood cells causedby Plasmodium falciparum. Blood, 74:855–861, 1989.112Bibliography[154] D. Needham and R. M. Hochmuth. Rapid flow of passive neutrophilsinto a 4 µm pipet and measurement of cytoplasmic viscosity. J.Biomech. Eng., 112:269–276, 1990.[155] H. Noguchi and G. Gompper. Fluid vesicles with viscous membranesin shear flow. Phys. Rev. Lett., 93:258102, 2004.[156] H. Noguchi and G. Gompper. Shape transitions of fluid vesicles andred blood cells in capillary flows. PNAS, 102:14159–14164, 2005.[157] I. L. Novak, B. M. Slepchenko, A. Mogilner, and L. M. Loew. Co-operativity between cell contractility and adhesion. Phys. Rev. Lett.,93:268109, 2004.[158] Y. Ohta, J. H. Hartwig, and T. P. Stossel. FilGAP, a Rho- and ROCK-regulated GAP for Rac binds filamin A to control actin remodelling.Nat. Cell Biol., 8:803–814, 2006.[159] B. Olofsson. Rho guanine dissociation inhibitors: pivotal molecules incellular signaling. Cell. Signal., 11:545–554, 1999.[160] T. Omori, T. Ishikawa, D. Barthe`s-Biesel, A.-V. Salsac, J. Walter,Y. Imai, and T. Yamaguchi. Comparison between spring networkmodels and continuum constitutive laws: Application to the large de-formation of a capsule in shear flow. Phys. Rev. E, 83:041918, 2011.[161] Y. K. Park, M. Diez-Silva, G. Popescu, G. Lykotrafitis, W. Choi, M. S.Feld, and S. Suresh. Refractive index maps and membrane dynamics ofhuman red blood cells parasitized by Plasmodium falciparum. PNAS,105:13730–13735, 2008.[162] C. Pasternak, J. A. Spudich, and E. L. Elson. Capping of surfacereceptors and concomitant cortical tension are generated by conven-tional myosin. Nature, 341:549–551, 1989.[163] C. S. Peskin. The immersed boundary methods. Acta Numerica,11:479–517, 2002.[164] A. P. Pirentis, E. Peruski, A. L. Iordan, and D. Stamenovic´. A modelfor stress fiber realignment caused by cytoskeletal fluidization duringcyclic stretching. Cell. Mol. Bioeng., 4:67–80, 2011.[165] I. V. Pivkin and G. E. Karniadakis. Accurate coarse-grained modelingof red blood cells. Phys. Rev. Lett., 101:118105, 2008.113Bibliography[166] M. Postma, L. Bosgraaf, H. M. Loovers, and P. J. M. Van-Haastert.Chemotaxis: signaling modules join hands at front and tail. CMBO,5:35–40, 2004.[167] J. Pourati, A. Maniotis, D. Spiegel, J. L. Schaffer, J. P. Butler, J. J.Fredberg, D. E. Ingber, D. Stamenovic, and N. Wang. Is cytoskeletaltension a major determinant of cell deformability in adherent endothe-lial cells. Am. J. Physiol., 274:C1283–C1289, 1998.[168] C. Pozrikidis. Axisymmetric motion of a file of red blood cells throughcapillaries. Phys. Fluids, 17:031503, 2005.[169] M. Prass, K. Jacobson, A. Mogilner, and M. Radmacher. Direct mea-surement of the lamellipodial protrusive force in a migrating cell. J.Cell Biol., 174:767–772, 2006.[170] P. Preira, M.-P. Valignat, J. Bico, and O. The´odoly. Single cell rheom-etry with a microfluidic constriction: Quantitative control of fric-tion and fluid leaks between cell and channel walls. Biomicrofluidics,7:024111, 2013.[171] A. Prosperetti and G. Tryggvason. Computational Methods for Mul-tiphase Flow. Cambridge University Press, Cambridge, 2007.[172] J. Qian, H. Liu, Y. Lin, W. Chen, and H. Gao. A mechanochemicalmodel of cell reorientation on substrates under cyclic stretch. PLoSOne, 8:e65864, 2013.[173] A. Rabodzey, P. Alcaide, F. W. Luscinskas, and B. Ladoux. Me-chanical forces induced by the transendothelial migration of humanneutrophils. Biophys. J., 95:1428–1438, 2008.[174] S. Ramanujan and C. Pozrikidis. Deformation of liquid capsules en-closed by elastic membranes in simple shear flow: large deformationsand the effect of fluid viscosities. J. Fluid. Mech., 361:117–143, 1998.[175] A. C. Reymann, R. Boujemaa-Paterski, J. L. Martiel, C. Gue´rin,W. Cao, H. F. Chin, E. M. De La Cruz, M. The´ry, and L. Blanchoin.Actin network architecture can determine myosin motor activity. Sci-ence, 336:1310–1314, 2012.[176] B. D. Rogers, R. A. Dalrymple, M. Go´mez-Gesteira, and A. J. C.Crespo. User guide for the parallelSPHysics code using MPI v2.0., 2011.114Bibliography[177] S. S. Rosenfeld, J. Xing, L. Chen, and H. Sweeney. Myosin IIB is un-conventionally conventional. J. Biol. Chem., 278:27449–27455, 2003.[178] R. J. Russell, S. L. Xia, R. B. Dickinson, and T. P. Lele. Sarcomeremechanics in capillary endothelial cells. Biophys. J., 97:1578–1585,2009.[179] Y. Sako, K. Hibino, T. Miyauchi, Y. Miyamoto, M. Ueda, andT. Yanagida. Single-molecule imaging of signaling molecules in liv-ing cells. Single Mol., 1:159–163, 2000.[180] J. W. Sanger, J. M. Sanger, and B. M. Jockusch. Differences in thestress fibers between fibroblasts and epithelial cells. J. Cell Biol.,96:961–969, 1983.[181] S. Sanyal, S. Ege´e, G. Bouyer, S. Perrot, I. Safeukui, E. Bischoff,P. Buffet, K. W. Deitsch, O. Mercereau-Puijalon, P. H. David, T. J.Templeton, and C. Lavazec. Plasmodium falciparum STEVOR pro-teins impact erythrocyte mechanical properties. Blood, 119:e1–e8,2011.[182] K. Sato, T. Adachi, M. Matsuo, and Y. Tomita. Quantitative evalua-tion of threshold fiber strain that induces reorganization of cytoskeletalactin fiber structure in osteoblastic cells. J. Biomech., 38:1895–1901,2005.[183] K. Sato, T. Adachi, M. Matsuo, and Y. Tomita. Measurement oflocal strain on cell membrane at initiation point of calcium signal-ing response to applied mechanical stimulus in osteoblastic cells. J.Biomech., 40:1246–1255, 2007.[184] M. Sato, D. R. Theret, L. T. Wheeler, N. Ohshima, and R. M.Nerem. Application of the micropipette technique to the measure-ment of cultured porcine aortic endothelial cell viscoelastic properties.J. Biomech. Eng., 112:263–268, 1990.[185] B. Schnitzer, T. Sodeman, M. L. Mead, and P. G. Contacos. Pit-ting function of the spleen in malaria: ultrastructural observations.Science, 177:175–177, 1972.[186] B. Schnitzer, T. M. Sodeman, M. L. Mead, and P. G. Contacos. Anultrastructural study of the red pulp of the spleen in malaria. Blood,41:207–218, 1973.115Bibliography[187] T. W. Secomb and J. F. Gross. Flow of red blood cells in narrowcapillaries: role of membrane tension. Int. J. Microcirc. Clin. Exp.,2:229–240, 1983.[188] T. W. Secomb and R. Skalak. A two-dimensional model for capillaryflow of an asymmetric cell. Microvasc. Res., 24:194–203, 1982.[189] T. W. Secomb, R. Skalak, N. O¨zkaya, and J. F. Gross. Flow of axisym-metric red blood cells in narrow capillaries. J. Fluid Mech., 163:405–423, 1986.[190] U. Seifert. Configurations of fluid membranes and vesicles. Adv. Phys.,46:13–137, 1997.[191] C. Y. Seow. Hill’s equation of muscle performance and its hiddeninsight on molecular mechanisms. J. Gen. Physiol., 142:561–573, 2013.[192] Y. M. Serebrennikova, J. Patel, W. K. Milhous, and L. H. Garc´ıa-Rubio. Quantitative analysis of morphological alterations in Plasmod-ium falciparum infected red blood cells through theoretical interpre-tation of spectral measurements. J. Theor. Biol., 265:493–500, 2010.[193] M. S. Shadloo, A. Zainali, M. Yildiz, and A. Suleman. A robust weaklycompressible SPH method and its comparison with an incompressibleSPH. Int. J. Numer. Methods Eng., 89:939–956, 2012.[194] S. Shao and E. Y. M. Lo. Incompressible SPH method for simulatingNewtonian and non-Newtonian flows with a free surface. Adv. WaterResour., 26:787–800, 2003.[195] J. P. Shelby, J. White, K. Ganesan, P. K. Rathod, and D. T. Chiu. Amicrofluidic model for single-cell capillary obstruction by Plasmodiumfalciparum-infected erythrocytes. PNAS, 100:14618–14622, 2003.[196] H. Shi, A. Li, J. Yin, K. S. W. Tan, and C. T Lim. AFM studyof the cytoskeletal structures of malaria infected erythrocytes. Proc.ICBME, 23:1965–1968, 2009.[197] Y. Shifrin, P. D. Arora, Y. Ohta, D. A. Calderwood, and C. A. Mc-Culloch. The role of FilGAP-Filamin A interactions in mechanopro-tection. Mol. Biol. Cell, 20:1269–1279, 2009.[198] J. W. Shin, A. Buxboim, K. R. Spinler, J. Swift, D. A. Christian, C. A.Hunter, C. Le´on, C. Gachet, P. C. Dave, P. Dingal, I. L. Ivanovska,116BibliographyF. Rehfeldt, J. A. Chasis, and D. E. Discher. Contractile forces sustainand polarize hematopoiesis from stem and progenitor cells. Cell StemCell, 14:81–93, 2014.[199] R. Skalak, A. Tozeren, R. P. Zarda, and S. Chien. Strain energyfunction of red blood cell membranes. Biophys. J., 13:245–264, 1973.[200] M. Soares e Silva, M. Depken, B. Stuhrmann, M. Korsten, F. C.MacKintosh, and G. H. Koenderink. Active multistage coarseningof actin networks driven by myosin motors. Proc. Natl. Acad. Sci.USA, 108:9408, 2011.[201] J. A. Spudich. The myosin swinging cross-bridge model. Nat. Rev.Mol. Cell Bio., 2:387–392, 2001.[202] M. R. Stachowiak and B. O’Shaughnessy. Recoil after severing revealsstress fiber contraction mechanisms. Biophys. J., 97:462–471, 2009.[203] Y. Sui, Y. T. Chew, P. Roy, X. B. Chen, and H. T. Low. Transient de-formation of elastic capsules in shear flow: Effect of membrane bendingstiffness. Phys. Rev. E, 75:066301, 2007.[204] Y. Sui, Y. T. Chew, P. Roy, Y. P. Cheng, and H. T. Low. Dynamicmotion of red blood cells in simple shear flow. Phys. Fluids, 20:112106,2008.[205] S. Suresh. Mechanical response of human red blood cells in health anddisease: Some structure-property-function relationships. J. Mater.Res., 21:1871–1877, 2006.[206] S. Suresh. Biomechanics and biophysics of cancer cells. Acta Mater.,55:3989–4014, 2007.[207] S.-Y. Tee, J. Fu, C. S. Chen, and P. A. Janmey. Cell shape andsubstrate rigidity both regulate cell stiffness. Biophys. J., 100:L25–L27, 2011.[208] D. P. Theret, M. J. Levesque, M. Sato, R. M. Nerem, and L. T.Wheeler. The application of a homogeneous half-space model in theanalysis of endothelial cell micropipette measurements. J. Biomech.Eng., 110:190–199, 1988.[209] O. Thoumine and A. Ott. Time scale dependent viscoelastic and con-tractile regimes in fibroblasts probed by microplate manipulations. J.Cell Sci., 110:2109–2116, 1997.117Bibliography[210] L. Tilley, M. W. A. Dixon, and K. Kirk. The Plasmodium falciparum-infected red blood cell. Int. J. Biochem. Cell Biol., 43:839–842, 2011.[211] H. P. Ting-Beall, D. Needham, and R. M. Hochmuth. Volume andosmotic properties of human neutrophils. Blood, 81:2774–2780, 1993.[212] A. Tondon, H.-J. Hsu, and R. Kaunas. Dependence of cyclic stretch-induced stress fiber reorientation on stretch waveform. J. Biomech.,45:728–735, 2012.[213] X. Trepat, L. Deng, S. S. An, D. Navajas, D. J. Tschumperlin, W. T.Gerthoffer, J. P. Butler, and J. J. Fredberg. Universal physical re-sponses to stretch in the living cell. Nature, 447:592–595, 2007.[214] M. A. Tsai, R. S. Frank, and R. E. Waugh. Passive mechanical behav-ior of human neutrophils: Power-law fluid. Biophys. J., 65:2078–2088,1993.[215] K. Tsubota, S. Wada, and T Yamaguchi. Particle method for computersimulation of red blood cell motion in blood flow. Comput. Meth.Programs Biomed., 83:139–146, 2006.[216] K. I. Tsubota and S. Wada. Elastic force of red blood cell membraneduring tank-treading motion: Consideration of the membrane’s natu-ral state. Inter. J. Mech. Sci., 52:356–364, 2010.[217] S. A. Vanapalli, M. H. G. Duits, and F. Mugele. Microfluidics as afunctional tool for cell mechanics. Biomicrofluidics, 3:012006, 2009.[218] B. Vanderlei, J. J. Feng, and L. Edelstein-Keshet. A computationalmodel of cell polarization and motility coupling mechanics and bio-chemistry. Multiscale Model. Simul., 9:1420–1443, 2011.[219] D. Vavylonis, J. Wu, S. Hao, B. O’Shaughnessy, and T. D. Pollard.Assembly mechanism of the contractile ring for cytokinesis by fissionyeast. Science, 319:97–100, 2008.[220] C. Veigel, J. E. Molloy, S. Schmitz, and K. Jones. Load-dependentkinetics of force production by smooth muscle myosin measured withoptical tweezers. Nat. Cell Biol., 5:980–986, 2003.[221] C. Veigel, S. Schmitz, F. Wang, and J. R. Sellers. Load-dependentkinetics of myosin-V can explain its high processivity. Nat. Cell Biol.,7:861–869, 2005.118Bibliography[222] M. Vicente-Manzanares, M. A. Koach, L. Whitmore, M. L. Lamers,and A. F. Horwitz. Segregation and activation of myosin IIB createsa rear in migrating cells. J. Cell Biol., 183:543–554, 2008.[223] M. Vicente-Manzanares, J. Zareno, L. Whilmore, C. K. Choi, andA. F. Horwitz. Regulation of protrusion, adhesion dynamics, andpolarity by myosins IIA and IIB in migrating cells. J. Cell Biol.,176:573–580, 2007.[224] S. Wada and R. Kobayashi. Numerical simulation of various shapechanges of a swollen red blood cell by decrease of its volume. Trans.Jpn. Soc. Mech. Eng., 69:14–21, 2003.[225] S. Walcott and S. X. Sun. A mechanical model of actin stress fiberformation and substrate elasticity sensing in adherent cells. PNAS,107:7757–7762, 2010.[226] F. Wang, M. Kova´cs, A. Hu, J. Limouze, E. V. Harvey, and J. R.Sellers. Kinetic mechanism of non-muscle myosin IIB: Functionaladaptations for tension generation and maintenance. J. Biol. Chem.,278:27439–27448, 2003.[227] I. Wang, A. Politi, N. Tania, Y. Bai, M. Sanderson, and J. Sneyd. Amathematical model of airway and pulmonary arteriole smooth mus-cle. Biophys. J., 94:2053–2064, 2008.[228] N. Wang, I. M. Tolic´-Nørrelykke, J. Chen, S. M. Mijailovich, J. P.Butler, J. J. Fredberg, and D. Stamenovic´. Cell prestress. I. Stiffnessand prestress are closely associated in adherent contractile cells. Am.J. Physiol.-Cell Ph., 282:C606–C616, 2002.[229] Q. Wang, J. J. Feng, and L. M. Pismen. A cell-level biomechanicalmodel of drosophila dorsal closure. Biophys. J., 103:2265–2274, 2012.[230] D. M. Warshaw, J. M. Desrosiers, S. S. Work, and K. M. Trybus.Smooth muscle myosin cross-bridge interactions modulate actin fila-ment sliding velocity in vitro. J. Cell Biol., 111:453–463, 1990.[231] Z. Wei, V. S. Deshpande, R. M. McMeeking, and A. G. Evans. Analysisand interpretation of stress fiber organization in cells subject to cyclicstretch. J. Biomech. Eng., 130:031009, 2008.[232] S. J. Winder and K. R. Ayscoughl. Cell science at a glance: Actin-binding proteins. J. Cell Sci., 118:651–654, 2005.119Bibliography[233] L. Wolff, P. Fernandez, and K. Kroy. Resolving the stiffening-softeningparadox in cell mechanics. PLoS One, 7:e40063, 2012.[234] K. Wong, O. Pertz, K. Hahn, and H. Bourne. Neutrophil polarization:Spatiotemporal dynamics of RhoA activity support a self-organizingmechanism. PNAS, 103:3639–3644, 2006.[235] T. Wu and J. J. Feng. Simulation of malaria-infected red blood cellsin microfluidic channels: Passage and blockage. Biomicrofluidics,7:044115, 2013.[236] T. Xu, D. Vavylonis, F. Tsai, G. H. Koenderink, W. Nie, E. Yusuf,I. Lee, J. Wu, and X. Huang. SOAX: A software for quantification of3D biopolymer networks. Sci. Rep., 5:9081, 2015.[237] B. Yap and R. D. Kamm. Cytoskeletal remodeling and cellular activa-tion during deformation of neutrophils into narrow channels. J. Appl.Physiol., 99:2323–2330, 2005.[238] B. Yap and R. D. Kamm. Mechanical deformation of neutrophilsinto narrow channels induces pseudopod projection and changes inbiomechanical properties. J. Appl. Physiol., 98:1930–1939, 2005.[239] A. Yazdani and P. Bagchi. Three-dimensional numerical simulationof vesicle dynamics using a front-tracking method. Phys. Rev. E,85:056308, 2012.[240] A. Yeung and E. Evans. Cortical shell-liquid core model for pas-sive flow of liquid-like spherical cells into micropipettes. Biophys. J.,56:139–149, 1989.[241] Z. Yu, X. Shao, and A. Wachs. A ctitious domain method for partic-ulate ows with heat transfer. J. Comput. Phys., 217:424–452, 2006.[242] Z. Yu, A. Wachs, and Y. Peysson. Numerical simulation of particlesedimentation in shear-thinning uids with a ctitious domain method.J. Non-Newtonian Fluid Mech., 136:126–139, 2006.[243] Y. Zhang, C. Huang, S. Kim, M. Golkaram, M. W. A. Dixon, L. Tilley,J. Li, S. Zhang, and S. Suresh. Multiple stiffening effects of nanoscaleknobs on human red blood cells infected with plasmodium falciparummalaria parasite. PNAS, Early Edition, 2015.120Bibliography[244] D. V. Zhelev and R. M. Hochmuth. Mechanically stimulated cytoskele-ton rearrangement and cortical contraction in human neutrophils. Bio-phys. J., 68:2004–2014, 1995.[245] C. Zhou, P. Yue, and J. J. Feng. Deformation of a compound dropthrough a contraction in a pressure-driven pipe flow. Int. J. MultiphaseFlow, 34:102 – 109, 2008.121


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items