ABSTRACT
Title of dissertation: MECHANICAL CHARACTERIZATION OF
NORMAL AND CANCEROUS BREAST
TISSUE SPECIMENS USING ATOMIC
FORCE MICROSCOPY
Rajarshi Roy
Dissertation directed by: Professor Jaydev P. Desai
Department of Mechanical Engineering,
Robotics, Automation, and Medical
Systems (RAMS) Laboratory.
Breast cancer is one of the most common malignancies among women world-
wide. Conventional breast cancer diagnostic methods involve needle-core biopsy
procedures, followed by careful histopathological inspection of the tissue specimen
by a pathologist to identify the presence of cancerous lesions. However, such in-
spections are primarily qualitative and depend on the subjective impressions of ob-
servers. The goal of this research is to develop approaches for obtaining quantitative
mechanical signatures that can accurately characterize malignancy in pathological
breast tissue. The hypothesis of this research is that by using contact-mode Atomic
Force Microscopy (AFM), it is possible to obtain differentiable measures of stiffness
of normal and cancerous tissue specimens.
This dissertation summarizes research carried out in addressing key experi-
mental and computational challenges in performing mechanical characterization on
breast tissue. Firstly, breast tissue specimens studied were 600 µm in diameter,
about six times larger than the range of travel of conventional AFM X-Y stages
used for imaging applications. To scan tissue properties across large ranges, a semi-
automated image-guided positioning system was developed that can be used to per-
form AFM probe-tissue alignment across distances greater than 100 µm at multiple
magnifications. Initial tissue characterization results indicate that epithelial tissue
in cancer specimens display increased deformability compared to epithelial tissue
in normal specimens. Additionally, it was also observed that the tissue response
depends on the patient from whom the specimens were acquired.
Another key challenge addressed in this dissertation is accurate data analysis
of raw AFM data for characterization purposes. Two sources of uncertainty typically
influence data analysis of AFM force curves: the AFM probe’s spring constant and
the contact point of an AFM force curve. An error-in-variable based Bayesian
Changepoint algorithm was developed to quantify estimation errors in the tissue’s
elastic properties due to these two error sources. Next, a parametric finite element
modeling based approach was proposed in order to account for spatial heterogeneity
in the tissue response. By using an exponential hyperelastic material model, it
was shown that it is possible to obtain more accurate material properties of tissue
specimens as opposed to existing analytical contact models.
The experimental and computational strategies proposed in this dissertation
could have a significant impact on high-throughput quantitative studies of bioma-
terials, which could elucidate various disease mechanisms that are phenotyped by
their mechanical signatures.
MECHANICAL CHARACTERIZATION OF NORMAL AND
CANCEROUS BREAST TISSUE SPECIMENS USING ATOMIC
FORCE MICROSCOPY
by
Rajarshi Roy
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2014
Advisory Committee:
Professor Jaydev P. Desai, Chair
Professor Carol Keefer, Dean’s Representative
Professor Balakumar Balachandran
Professor Nikhil Chopra
Professor Teng Li
c© Copyright by
Rajarshi Roy
2014
1. Portions reprinted, with permission, from: Rajarshi Roy, Wenjin Chen,
Lauri Goodell, Jun Hu, David J. Foran and Jaydev P. Desai, “Microarray-
facilitated Mechanical Characterization of Breast Tissue Pathology Samples
using contact mode Atomic Force Microscopy (AFM)”, International Confer-
ence on Biomedical Robotics and Biomechatronics pp. 710 - 715, 2010.
2. Portions reprinted, with permission, from: Rajarshi Roy, Wenjin Chen, Lei
Cong, Lauri Goodell and David J. Foran, Jaydev P. Desai, “Towards Auto-
mated Mechanical Characterization of Biological Samples using Atomic Force
Microscopy (AFM)”, 2011 ASME Dynamic Systems and Control Conference,
November 2011.
3. Portions reprinted, with permission, from: Rajarshi Roy, Lei Li, Carol L.
Keefer and Jaydev P. Desai, “An Error-In-Variables (EIV) based Bayesian
Probabilistic Approach to Estimating Cell Mechanical Properties using Atomic
Force Microscopy”, International Conference on Biomedical Robotics and Biomecha-
tronics, pp. 389–394, June 2012.
4. Portions reprinted, with permission, from: Rajarshi Roy, Wenjin Chen, Lei
Cong, Lauri Goodell, David J. Foran and Jaydev P. Desai, “A Semi-Automated
Positioning System for contact-mode Atomic Force Microscopy (AFM)”, IEEE
Transactions on Automation Science and Engineering, vol. 10 no. 2, pp. 462-
465, 2013.
5. Portions reprinted, with permission, from: Rajarshi Roy, Wenjin Chen, Lei
Cong, Lauri Goodell, David J. Foran and Jaydev P. Desai, “Probabilistic
Estimation of Mechanical Properties of Biomaterials Using Atomic Force Mi-
croscopy”, IEEE Transactions on Biomedical Engineering, vol. 61 no. 2, pp.
547-556, 2014.
Acknowledgments
I am deeply grateful to all those people who made this dissertation possible
and made my graduate school experience something that I would cherish for the
rest of my life.
Firstly, I would like to thank my advisor, Professor Jaydev P. Desai for giving
me the opportunity to work in his laboratory on an extremely challenging project
over the past five years. His enthusiasm, constant encouragement and patience,
especially during my initial graduate school years, is something I am profoundly
grateful for. He always made himself available, even at late nights, for troubleshoot-
ing with the AFM. It has been a pleasure to work and learn from him over the past
five years.
I would also like to thank my collaborators, Dr. Wenjin Chen and Professor
David J. Foran from the Rutgers Cancer Institute of New Jersey. Without their
critical feedback at various stages of this project, this dissertation would not have
been possible. It has been a pleasure to work and collaborate with them.
Next, I would like to thank Dr. Carol Keefer taking the time out of her busy
schedule to teach me various intricacies of light microscopy and the components of
optical microscopes. I could always count on her to provide me with lab accessories
urgently whenever I ran out of them.
My colleagues at the Robotics, Automation and Medical Systems lab have
been instrumental in enriching my graduate experience - I learnt a lot from them.
Dr. Anand Pillarisetti helped me get acquainted with the AFM. Dr. U-Xuan Tan
ii
provided great help when I was starting with the MP-285 micromanipulator. Dr.
Mingyen Ho greatly assisted me with any computer hardware problems I had. I had
a lot of fruitful discussions on computational modeling with Dr. Kevin Lister. It
has been a pleasure discussing various aspects of this project with Dr. Yang Bo,
Elif Ayvali, Chad Kessens, Hyuntae Kim and Dr. Hardik Pandya. I would like to
express my gratitude to all members of the RAMS laboratory who have helped me
during my graduate studies.
Special thanks to Mr. Melvin Fields for his prompt response with any com-
puter related problems I had.
I would also like to thank Xavier Pelle, Suvajyoti Guha, Snehaunshu Chowd-
hury, Ishita Chakraborty, Sandip Biswas and Sagar Chowdhury for their friendship
during these last five years.
I owe my deepest gratitude to my parents, who have constantly stood beside
me through all the difficult phases I had during these last five years. The values and
morals that I have imbibed from them kept me strong when the going got tough.
Thank you!
I would also like to acknowledge the financial support from the National Sci-
ence Foundation (NSF) and National Institutes of Health (NIH) for the research
presented in this dissertation.
Lastly, I would like to apologize to those I’ve inadvertently left out. Thanks
again!!
iii
Table of Contents
List of Tables vi
List of Figures vi
List of Abbreviations xii
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Proposed Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Pilot Mechanical Characterization Studies on Breast Tissue Pathology sam-
ples using Atomic Force Microscopy 10
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Material and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Tissue Microarray (TMA) preparation and Annotation . . . . 10
2.2.2 AFM experimental setup . . . . . . . . . . . . . . . . . . . . . 14
2.2.3 Contact-mode operation of AFM . . . . . . . . . . . . . . . . 16
2.2.4 Operating modes for AFM indentation . . . . . . . . . . . . . 19
2.2.4.1 Deflection-controlled AFM force curves . . . . . . . . 19
2.2.4.2 Indentation-controlled AFM force curves . . . . . . . 21
2.2.5 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3 Improvements to the AFM experimental protocol for tissue characterization 33
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2 Problem Statement and Proposed Solution . . . . . . . . . . . . . . . 34
3.3 AFM experimental setup with the IGPS . . . . . . . . . . . . . . . . 39
3.4 Image-Guided Navigation . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4.1 Tracking of the ROI . . . . . . . . . . . . . . . . . . . . . . . 40
3.4.2 Control Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.5 ROI positioning protocol . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.6 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6.1 Characterization of Epithelium-Stroma boundary . . . . . . . 49
3.6.2 Characterization of patient specific tissue specimens . . . . . . 53
3.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4 Probabilistic Estimation of the Elastic Modulus 58
4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2 Problem Statement and Proposed Solution . . . . . . . . . . . . . . . 59
4.2.1 Uncertainty in the contact Point . . . . . . . . . . . . . . . . . 59
4.2.2 Uncertainty in AFM Probe Calibration . . . . . . . . . . . . . 64
iv
4.3 Bayesian Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.1 Transformation of raw data . . . . . . . . . . . . . . . . . . . 71
4.3.2 Classical Error-In-Variables Model . . . . . . . . . . . . . . . 75
4.3.3 EIV-based Changepoint Modeling: Likelihood . . . . . . . . . 77
4.3.4 EIV-based Changepoint Modeling: Posterior . . . . . . . . . . 81
4.3.5 EIV-based Changepoint Modeling: Gibbs Sampling . . . . . . 82
4.4 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4.1 Implementation of the Gibbs Sampler . . . . . . . . . . . . . . 87
4.4.2 Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.4.3 AFM studies on histological breast tissue . . . . . . . . . . . . 92
4.4.4 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . 96
4.4.4.1 Varying σsp . . . . . . . . . . . . . . . . . . . . . . . 96
4.4.4.2 Varying µsp . . . . . . . . . . . . . . . . . . . . . . . 98
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5 Constitutive material modeling of the tissue specimens under AFM indenta-
tion 104
5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.2 Experimental Challenges in enforcing Hertzian assumptions . . . . . . 105
5.3 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.4 Problem Statement and Proposed Solution . . . . . . . . . . . . . . . 116
5.5 Finite Element Modeling . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.5.1 Discretization and Parameter Estimation . . . . . . . . . . . . 121
5.5.2 Constitutive Material Models . . . . . . . . . . . . . . . . . . 125
5.5.3 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . 129
5.6 Contact Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
5.6.1 BDP approach without continuity [1] . . . . . . . . . . . . . . 134
5.6.2 Continuity-constrained weighted BDP approach . . . . . . . . 136
5.6.3 Choice of fi for automated data analysis . . . . . . . . . . . . 139
5.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
5.7.1 Parameter reidentification on simulated force curves . . . . . . 141
5.7.2 Implementation on elastic map . . . . . . . . . . . . . . . . . 145
5.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
6 Conclusions and Future Work 151
6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
6.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
6.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
A Finite indentation response of Hertz and Dimitriadis contact models 156
Bibliography 158
v
List of Tables
2.1 AFM probe parameters used for characterization experiments . . . . 16
2.2 AFM probe spring constant for all experiments . . . . . . . . . . . . 26
2.3 Summary of AFM indentation experiments. Results shown as mean
± standard deviation. P-values less than α(= 0.05) are highlighted
in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.1 X-Y travel range of automated AFM scanners, arranged in increasing
order of their travel range. . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2 AFM experimental setting parameters. . . . . . . . . . . . . . . . . . 53
3.3 Summary of AFM raster-scanning results on patient specific tissue
specimens. Results expressed as mean ± standard deviation. . . . . . 54
4.1 Variation of elastic modulus with change in contact point. . . . . . . 61
4.2 Summary of AFM probe calibration results using the thermal method.
Results expressed as mean ± standard deviation. . . . . . . . . . . . 69
4.3 Conjugate Prior Distributions . . . . . . . . . . . . . . . . . . . . . . 81
4.4 Hyperparameter Values . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.5 Model performance on simulated data . . . . . . . . . . . . . . . . . . 92
5.1 Summary of implementation of the EIV-Bayesian Changepoint algo-
rithm on two experimental force curves E1 and E2 (Figs. 5.5 and 5.6
respectively acquired on a tissue specimen.) . . . . . . . . . . . . . . 117
5.2 Parameters used to construct the simulated force curves. . . . . . . . 141
5.3 Summary of algorithm performance on force curves S1-S3. Units of
E and E0 are in kPa. G.T. indicates ground truth. . . . . . . . . . . 145
List of Figures
1.1 Stained histopathological images of (a) ductal carcinoma in situ, (b)
lobular carcinoma in situ and (c) invasive ductal carcinoma [2]. . . . . 2
1.2 Schematics of (a) Micropipette Aspiration, (b) Magnetic Twisting
Cytometry, (c) Optical Tweezer and (d) contact-mode Atomic Force
Microscopy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Length scales of various biological specimens, together with the ap-
propriate force measuring instruments. . . . . . . . . . . . . . . . . . 6
2.1 (a) Stained slide image and (b) adjacent unstained slide image c©2010
IEEE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
vi
2.2 An example of one tissue sample as has been used in different steps of
experiment: (a) A low resolution stained image of the tissue microar-
ray. (b) The VM image of an H&E stained tissue core. (c) The same
core annotated by a pathologist. (d) The unstained adjacent tissue
core (VM image). Please note that tissue distribution in this image
is nearly the same as (b) but nearly impossible to distinguish by eye.
(e) Image taken from the AFM microscope during AFM experiments
c©2013 IEEE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 AFM experimental setup c©2010 IEEE. . . . . . . . . . . . . . . . . . 14
2.4 Hydrated tissue specimen on microscope slide. . . . . . . . . . . . . . 15
2.5 CAD drawing of a typical AFM probe used for tissue characterization
experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.6 Schematic of contact-mode AFM operation. . . . . . . . . . . . . . . 17
2.7 (a) Screenshot of MFP3D software panel used to acquire AFM force
curves in the deflection-controlled mode. (b) Acquired AFM force
curve with a deflection setpoint of 50 nm. The red and blue arrows
indicate the direction of the probe’s motion. . . . . . . . . . . . . . . 20
2.8 Screenshot of cropped MFP3D software panel used to acquire AFM
force curves in the indentation-controlled mode. (a) Screenshot of 10
nm deflection-trigger. (b) Acquired force curve with 10 nm deflection
trigger. (c) Screenshot of 250 nm indentation trigger. (d) Acquired
force curve with 250 nm indentation trigger. The red and blue arrows
indicate the direction of the probe’s motion. . . . . . . . . . . . . . . 22
2.9 Contact estimation algorithm. . . . . . . . . . . . . . . . . . . . . . . 24
2.10 Hertzian contact of an elastic half-space by a sphere of radius R. . . . 25
2.11 Histology images (20X magnification) cropped from four representa-
tive regions probed in set 1. (a) and (b) are regions with normal
epithelial (A8) and cancerous epithelial (A13) tissue respectively, (c)
and (d) are regions with normal stromal (A22) and cancerous stromal
(A12) tissue respectively c©2010 IEEE. . . . . . . . . . . . . . . . . . 26
2.12 Force-Indentation plots obtained in set 1 at (a) normal epithelial re-
gions (A8), (b) cancer epithelial regions (A13), (c) normal stromal
regions (A22) and (d) cancer stromal regions (A12) along with their
Hertzian fits and the mean R2 values. Inside each annotated region,
roughly 5-8 force curves were taken. Note that despite a target in-
dentation of 250 nm, the tissue indentation varied from 150 - 500 nm.
This was primarily on account of incorrect contact estimation in the
indentation-controlled mode. . . . . . . . . . . . . . . . . . . . . . . . 27
3.1 Flowchart describing manual AFM indentation experiments. (M) = Man-
ual, (A) = Automated. . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
vii
3.2 Manual registration between unstained and stained images. (a) Stained
TMA core image annotated by the pathologist as regions A1, A2 and
A3. (b) Brightfield image while performing AFM indentation on a
point in annotation A3 after manually registering the stained image
(c) Stained TMA core image overlaid with the points probed during
AFM experiments on the unstained slide. Observe that annotation
A2 has been incorrectly sampled by the AFM. . . . . . . . . . . . . . 36
3.3 AFM Experimental Setup with the MP-285 micromanipulator c©2013
IEEE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4 Extracted features with strong intensity gradients. Observe that
prominent features are impurities floating in the liquid environment
of the sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.5 Normalized cross correlation algorithm used for tracking at two time
instants (a) and (b). The top images correspond to It = I(px, py, t)
with the ROI marked with a white box, while the bottom ones are the
result images R, where with the estimated ROI position (pˆx,t+1, pˆy,t+1)
appear as spots of high intensity. . . . . . . . . . . . . . . . . . . . . 42
3.6 Calibration of the scaling matrix S using a cover-slip with uniform
grids at 100 µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.7 Tracking protocol for probe-ROI alignment at 10X and 20X. (a)
Stained image of tissue core annotated by the pathologist at low
resolution (10X). (b) Brightfield image from the AFM-optical mi-
croscope. A coarse ROI is selected to match one of the annotated
regions in (a). (c) Positioning at 10X. (d) Registration of the coarse
ROI at 20X c©2013 IEEE. . . . . . . . . . . . . . . . . . . . . . . . . 45
3.8 Flowchart describing AFM indentation experiments using the image-
guided positioning system (IGPS). (M) = Manual, (A) = Automated. 47
3.9 (a) Representative AFM force curve on a tissue sample mounted on
the end-effector. (b) Corresponding force-indentation curve overlaid
with the Hertzian fit and (c) Tracking Performance c©2013 IEEE. . . 48
3.10 Elastic maps on normal tissue cores along with their correponding
brightfield images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.11 Elastic maps on cancer tissue cores along with their correponding
brightfield images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.12 Uninterpretatble force curve with no prominent non-contact regime. . 52
4.1 Representative AFM force curve on mESC, with candidate contact
points shown in red blocks: (a) Whole AFM force curve and (b)
transition region from non-contact to contact regime. . . . . . . . . . 60
4.2 AFM Probe Calibration using the thermal method. (a) An AFM force
curve on a tissue-free hard glass surface and (b) thermal spectrum of
the probe in PBS solution. The rectangular probe of Table 4.2 (rated
compliance of 0.222 m/N) was used to generate these results c©2014
IEEE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
viii
4.3 Representative AFM force curve of d v/s z and d v/s δ. Please note
that the d v/s δ curve is displaced 10nm on the deflection axis for
ease of visualization c©2014 IEEE. . . . . . . . . . . . . . . . . . . . . 72
4.4 Schematic of the EIV-Bayesian Changepoint algorithm c©2014 IEEE. 82
4.5 Results of the EIV-Bayesian Changepoint algorithm on simulated
data with n = 1000 datapoints. (a) Simulated AFM force curve with
contact occurring at datapoint k = 650 and (b) marginal posterior of
the contact point index k, along with the posterior mean bkˆc = 653
and the 95% confidence intervals [642, 660]. . . . . . . . . . . . . . . . 89
4.6 Marginal posterior distributions in (a) β22 (b) s (c) σ1 and (d) σ2. . 91
4.7 Image of tissue region where probing was conducted to obtain the
results in this section. AFM indentation experiments on breast tissue
were conducted in the probing ROI c©2014 IEEE. . . . . . . . . . . . 93
4.8 Results of the EIV-Bayesian Changepoint algorithm on a representa-
tive AFM force curve on breast tissue with n = 783 datapoints. (a)
AFM force curve with the posterior mean of the contact point index
at bkˆc = 433. (b) Marginal posterior of k. . . . . . . . . . . . . . . . 95
4.9 (a) Marginal posterior in s and (b) marginal Posterior of E. . . . . . 96
4.10 Boxplots of the marginal posterior distributions (a) in E, (b) in s and
(c) k. Please note that the σsd has been assumed to be equal to σsp.
The fourth boxplot (σsp = 0.032) corresponds to the experimental
probe compliance variance as obtained in Table 4.2 c©2014 IEEE. . . 97
4.11 Boxplots indicating the effect of varying µsp on the marginal posteri-
ors (a) in E and (b) in s c©2014 IEEE. . . . . . . . . . . . . . . . . . 100
5.1 Schematic of force-displacement behavior for (a) adhesion-free inden-
tation, (b) indentation with adhesion and (c) representative AFM
force curve on breast tissue specimens showing absence of prominent
adhesive region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.2 Finite-indentation of a thin specimen by a rigid hemispherical indenter.107
5.3 Effect of heterogeneity in the spatial stiffness distribution. (a) An
elasticity map in the cancer core [from Fig. 3.11(b)], (b) correspond-
ing 2-D representation of the R-square of Hertzian fits and (c,d) repre-
sentative force-indentation plots showing the inadequacy of the Hertz
model to describe the entire ROI’s constituent behavior. . . . . . . . 110
5.4 Effect of finite deformation on the force response during spherical in-
dentation. (a) Simulated force response on a finite element model of a
tissue specimen (cyan) and a Hertzian force response of the same elas-
tic modulus (purple) and (b) corresponding indentation stress (σ∗) -
indentation strain (∗) using Eqns. (5.3a - c). . . . . . . . . . . . . . . 114
5.5 Results of Gibbs Sampling on the posterior distribution using (a)
Hertz contact model, (b) Dimitriadis’s contact model and (c) Long’s
contact model on force curve E1. . . . . . . . . . . . . . . . . . . . . 118
ix
5.6 Results of Gibbs Sampling on the posterior distribution using (a)
Hertz contact model, (b) Dimitriadis’s contact model and (c) Long’s
contact model on force curve E2. . . . . . . . . . . . . . . . . . . . . 119
5.7 Evidence of lack of material homogeneity in the direction of inden-
tation. (a) AFM force curves. Please note that the force curves
have been displaced to have the same approximate contact point at
(0,0).(b) Corresponding force-indentation curves. Note the change
in the nature of the force-indentation curves at 200 nm and 300 nm
indentation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.8 Finite element models: (a) undeformed mesh, (b) deformed mesh. . . 123
5.9 Flowchart of the algorithm used to obtain optimized material prop-
erties using the inverse FE approach. . . . . . . . . . . . . . . . . . . 124
5.10 Performance of the (a) Yeoh and (b) Fung hyperelastic models in
describing the force response for force curve shown in Fig. 5.6. . . . . 127
5.11 Effect of varying the hyperelastic parameters of the Yeoh (a,b,c) and
Fung (d,e) hyperelastic models individually. . . . . . . . . . . . . . . 129
5.12 Sensitivity of the Fung model parameters at varying indentations.
(a)∆ = 100 nm, (b) ∆ = 200 nm, (c) ∆ = 400 nm, (d) ∆ = 600 nm,
(e) ∆ = 900 nm and (f) ∆ = 1400 nm. Also shown in the third row
(g-i) are enlarged figures of (d-f) where the range of b is reduced to
1−6. The white dotted line at 400 kPa indicates the elastic modulus
of the linearly elastic material. . . . . . . . . . . . . . . . . . . . . . . 130
5.13 Classification curve, showing the variation of bmaxR2FE,thresh(∆) with ∆
for different values of R2FE,thresh. Also shown are the curve fits for
bmaxR2FE,thresh(∆). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
5.14 Performance of the three contact estimation methods for simulated
force curves S1 - S3 (shown row-wise). The first column shows the
sum of squared residuals (SSR) using unconstrained BDP (red) for
force curves S1-S3 (a),(d),(g) respectively. The second column shows
SSR variations using constrained BDP (blue and green) for S1-S3
(b),(e),(h) respectively. The third column shows the corresponding
contact point on the force curves S1-S3 (c),(f),(i) respectively. The
contact point estimates using the unconstrained and constrained BDP
approaches are also shown at the base of the SSR curve. . . . . . . . 143
5.15 Implementation of the proposed algorithm for automated characteri-
zation for simulated force curves S1 - S3 (shown row-wise). Note that
the inverse FE correction step was not required for curve S1. . . . . . 144
x
5.16 Implementation of the proposed algorithm for automated character-
ization on a tissue ROI shown previously in Fig. 5.3. (a) Elastic
map using the Dimitriadis model only without any Inverse FE cor-
rections, (b) Corresponding post-contact R2Dim fit, (c) Elastic map
(Eˆ) with inverse-FE correction, (d) Corresponding R2net fit, where
R2net = max(R2Dim, R2FE). White spaces indicate force curves that
were not analyzed due to presence of multiple potential contact points
such as those in Fig. 5.7. . . . . . . . . . . . . . . . . . . . . . . . . . 146
5.17 (a) Variation of estimated Fung parameter bˆ from the Inverse FE
analysis as a function of the estimated indentation (using the contact
estimate), (b) zoomed image. Also overlaid in red is bmax0.992(∆). . . . . 148
5.18 Representative force-indentation plot with low estimated b value (Eˆ0, bˆ)
= (153.21 kPa,0.005) with R2FE = 0.977. The Dimitriadis solution
was E¯Dim was 146.14 kPa with R2Dim = 0.975. . . . . . . . . . . . . . 148
A.1 (a) Variation of EHertz/E and EDim/E for R = 2.5 µm, (b) Corre-
sponding R2 fit quality, (c) Variation of EHertz/E and EDim/E for
R = 5.0 µm and (d) Corresponding R2 fit quality. . . . . . . . . . . . 157
xi
List of Abbreviations
AFM Atomic Force Microscopy
BDP Bi-domain Polynomial
CAD Computer-Aided Diagnostic
C.I. Confidence Interval
cov Covariance
DCIS Ductal carcinoma in situ
ECM Extracellular Matrix
EIV Error-In-Variables
FE Finite Element
FIEL Force Integration to Equal Limits
G.T. Ground Truth
H&E Hematoxylin and Eosin
IGPS Image-Guided Positioning System
InvOLS Inverse Optical Lever Sensitivity
i.i.d independent and identically distributed
KLT Kanade-Lucas-Tomasi
LCIS Lobular carcinoma in situ
LDV Laser Doppler Vibrometry
MA Micropipette Aspiration
MCMC Markov Chain Monte Carlo
mESC mouse Embryonic Stem Cell
MTC Magnetic Twisting Cytometry
OT Optical Tweezers
PBS Phosphate Buffered Saline
ROI Region-Of-Interest
SIFT Scale-Invariant Feature Transform
SSR Sum of Squared Residuals
TMA Tissue Microarray
thresh threshold
var Variance
VM Virtual Microscopy
NSF National Science Foundation
NIH National Institutes of Health
RAMS Robotics, Automation, and Medical Systems Laboratory
xii
Chapter 1
Introduction
1.1 Motivation
It is estimated by the American Cancer Society that 232,340 new cases of
invasive breast cancer are expected to be diagnosed among women in the US during
2013, and an estimated 40,030 breast cancer deaths alone are expected in 2013 [3].
Traditional assessments of cancer involve a biopsy procedure, followed by detailed
histopathological analysis of the sampled tissue to identify the presence of cancerous
lesions.
Histology is the branch of biology that studies the microscopic structure of
animal or plant tissues. The most common histological procedure is to embed fixed
tissue chunks in paraffin blocks, thinly section them, treat the sectioned tissue slices
with various staining techniques, and then study them under an optical microscope.
This method has enabled pathologists and clinicians to study human tissue in its
normal and diseased forms for hundreds of years and is still playing an indispensable
role in clinical pathology and biology research. Four types of tissue are present in
human body: epithelial tissue, stromal tissue, nervous tissue and muscles; and they
interweave together in distinctive patterns throughout the human body.
Most breast cancers are broadly classified by pathologists into the following
categories [2]:
1
1. Carcinoma in situ: These involve neoplastic proliferation that is limited to the
ducts [called ductal carcinoma in situ (DCIS) - see Fig. 1.1(a)] and lobules
[called lobular carcinoma in situ (LCIS) - see Fig. 1.1(b)] by the basement
membrane.
2. Invasive/Infiltrative carcinoma: These involve neoplastic proliferation that
have penetrated through the basement membrane into the stromal tissue [see
Fig. 1.1(c)].
(a) (b)
(c)
Figure 1.1: Stained histopathological images of (a) ductal carcinoma in situ, (b)
lobular carcinoma in situ and (c) invasive ductal carcinoma [2].
While the interpretation of stained images leads to insights into the onset and
progression of malignancy in a tissue specimen, this exercise is largely qualitative
and based on the subjective impressions of observers. To address this, researchers
have proposed a range of computer-aided diagnostic (CAD) methods that assist
physicians in rendering more informed clinical decisions. A comprehensive review of
2
the recent advances in automated interpretation of digital microscopy images and
CAD methods with relevance to histopathological tissue is given in [4].
The correlation between pathophysiological processes and the mechanical prop-
erties of the tissue they affect has been the subject of substantial research for a
number of years. Pathophysiological changes that occur in macroscale tissue have
been studied by imaging tools like ultrasound [5], magnetic resonance imaging [6]
and computed tomography [7], and there is reason to believe that these changes may
manifest at the single cellular level [8]. Li et al. [9] reported increased deformabil-
ity in cancerous breast epithelial cells compared to normal breast epithelial cells,
while similar increase in deformability was reported in cancerous bladder cells [10].
It has been hypothesized that the cellular cytoskeleton plays a dominant role in
the onset and progression of cancer. Actin reorganization was speculated to be the
cause of alteration in mechanical properties in bladder cells [10], breast epithelial
cells [9], mouse fibroblast cells [11] and keratinocytes [12]. Changes in deformability
in the extracellular matrix (ECM), that provides structural support to the cells,
was reported to activate malignancy in breast epithelial cells [13], while the role of
cell-ECM binding in cancer has been investigated in [14].
While the structural changes in cellular architecture and the surrounding ECM
induced by cancer has been widely documented at the cellular scale [8], the mor-
phological changes inside a developing tumor lesion is not clearly understood yet.
Indeed, the difference between cellular and tissue-level manifestation of cancer can
be noted by observing the apparent dichotomous impact of cancer across length
scales: pathologists use breast palpation to locate the presence of stiff nodules,
3
which correspond to hardening of the neighboring tumor stromal tissue [15], while
cell biophysicists widely report increased deformability in malignant cells compared
to healthy ones, as mentioned earlier [8].
There are several reasons for investigating the mechanical changes that accom-
pany onset and progression of cancer at the tissue level as opposed to the cellular
level. Studying the tumor at the tissue level provides direct insight into the underly-
ing architectural changes that occur within a developing lesion and the surrounding
tissues, and the degree of malignancy can be quantified in the tissue’s native envi-
ronment. In contrast, the relevance of single-cell study is somewhat questionable
since it is not representative of a three-dimensional tissue environment.
Another distinct advantage of studying tumor-level breast cancer is the im-
proved throughput associated with histological tissue preparation. Using techniques
like tissue microarray (TMA) [16] technology, where representative samples from
numerous biopsy procedures can be assembled onto a single microscope slide, it
is possible to characterize many sub-classes of cancer in a reasonably short period
of time. It is therefore enticing to imagine a battery of mechanical signatures of
various subclasses of cancer complementing the digitized stained-image information
from an imaging system to greatly enhance the capability of current computer-aided
diagnostic methods in making informed clinical decisions.
Advances in cancer biomechanics research has been supplemented by a surge in
the development of mechanical property measurement techniques at the micro/nano
scale [17]. Some of the most common methods used to quantify mechanical proper-
ties of biomaterials are given as follows (see Fig. 1.2 [17]):
4
1. Micropipette Aspiration (MA): In this technique, a biological cell is drawn into
a pipette using a known suction pressure. The aspiration length quantifies the
deformability of various types of cells.
2. Magnetic Twisting Cytometry (MTC): This method involves attaching mag-
netic beads to functionalized surfaces of cells. Application of a known magnetic
field deforms the surface of the cell, which characterizes the cellular elastic and
viscous properties.
3. Optical Tweezers (OT): In this method, a highly focused laser beam is aimed
at a dielectric bead, usually attached to the biological specimen. The high
intensity of the laser creates a “trap” which holds the bead and the attached
specimen at the focusing point. In many cases, two beads are attached to two
diametrically opposite surfaces of a cell.
4. Atomic Force Microscopy (AFM): The AFM system consists of a micro-cantilever
probe that is piezoelectrically controlled to deform tissue samples by a fixed
amount. Based on the AFM probe and tissue-sample interaction forces, the
cantilever deflects, which is then optically sensed. This deflection and the cor-
responding change in light intensity is related to the stiffness of the sample
probed by the AFM.
While the MA, MTC and OT based mechanical assays all possess the ability
to quantify mechanical properties of micro/nano scale biomaterials, AFM surpasses
these techniques in its versatility in quantifying mechanical properties of biomateri-
5
(a) (b)
(c)
Suction
Laser
Bead Bead
Tether
Bead
X-Y Stage
Z-Piezo
Laser SourcePhotodiode
(d)
Figure 1.2: Schematics of (a) Micropipette Aspiration, (b) Magnetic Twisting Cy-
tometry, (c) Optical Tweezer and (d) contact-mode Atomic Force Microscopy.
cells
tissue organ/tissue
proteins
Length [m]
DNA base
pair length
AFM
OT
MA
MTC
Figure 1.3: Length scales of various biological specimens, together with the appro-
priate force measuring instruments.
6
als across various length scales (see Fig 1.3). It also allows minimal sample prepara-
tion and therefore samples can be studied in their physiological environment. Apart
from its ability to measure nanoscale to picoscale forces, it is also widely used as an
imaging tool [18].
With this background in breast cancer and the mechanical assays that hold
the potential to quantify it, it is hypothesized that using the AFM in contact mode,
differentiable measures of local stiffness in histopathological breast tissue could re-
veal insights into the nature and progression of breast cancer that could supplement
traditional staining based histological assessments for diagnostic purposes.
1.2 Proposed Solution
The research objectives that have been addressed in this work are stated as
follows:
1. Tissue preparation protocol and preliminary feasibility study: In order to en-
sure that AFM sampling in the tissue is restricted to a region of biological
interest, it is necessary to have a predefined map of the tissue specimen which
can be used as a reference to navigate through the specimen under AFM inden-
tation. Virtual-microscopy (VM) enabled TMA technology for this purpose.
The results of preliminary characterization studies show that cancerous ep-
ithelial tissue exhibits increased deformability compared to normal epithelial
tissue.
2. Experimental improvements in the AFM sampling procedure: Human breast
7
tissue being highly heterogeneous wherein functional breast epithelial cells
are surrounded by stromal tissue in a complex intertwined manner, sampling
a particular cell type is highly labor-intensive and inefficient. To enhance
the positioning throughput prior to AFM indentation, a long-range image-
guided positioning system (IGPS) has been developed that improves efficiency
in AFM tissue characterization experiments significantly, and the developed
setup serves as a useful supplement to commercial AFM X-Y stages with
limited scanning range.
3. Error quantification in tissue stiffness estimates from indentation data: A crit-
ical pre-requisite to effective tissue characterization is accurate data analysis
of raw AFM data. Variability in two key experimental parameters give rise
to inaccurate estimation of the tissue mechanical properties: (1) the contact-
point in an AFM force curve, which relates to the probe’s vertical position
when it first makes contact with the tissue and (b) the probe’s spring con-
stant, which converts the probe deflection into the contact force. To this end,
an error-in-variables (EIV) based Bayesian Changepoint algorithm has been
developed to investigate the impact of these uncertainties in the mechanical
characterization process.
4. Constitutive modeling of breast tissue: The direct implication of spatial het-
erogeneity in breast tissue specimens is that tissue response might not be
adequately captured through linearized analytical contact models such as the
Hertz contact model. Violations to Hertzian contact theory has been investi-
8
gated by using numerical procedures and an exponential type phenomenolog-
ical hyperelastic constitutive model has been proposed to estimate nonlinear
tissue properties from AFM data.
1.3 Organization
Single-cell mechanical characterization using AFM has been well-researched
in the past[19] and this dissertation commences by following a cell-based AFM
characterization protocol as a starting point. In chapter 1, the background and
motivation for this study is discussed. In chapter 2, the operation of the AFM, the
tissue preparation and AFM experimental protocol are introduced and the results of
preliminary AFM indentation experiments on breast tissue specimens are discussed.
Following this, the development of a semi-automated positioning system is reported
in chapter 3, which improves efficiency in the AFM characterization experiments
on tissue compared to preliminary tissue characterization. In chapter 4, details
of a probabilistic mathematical model using EIV-Bayesian Changepoint scheme is
elaborated, which is used to estimate tissue mechanical properties from AFM data
after taking into consideration the innate errors in the AFM system. In chapter 5,
constitutive material modeling of the tissue specimens under AFM indentation is
discussed. Finally, in chapter 6, contributions of this dissertation and the future
goals of this project are discussed.
9
Chapter 2
Pilot Mechanical Characterization Studies on Breast Tissue
Pathology samples using Atomic Force Microscopy
2.1 Overview
This chapter discusses the results of pilot studies of AFM indentation on hu-
man breast tissue. Epithelial and stromal tissue regions were selected from normal
and cancerous tissue specimens and AFM indentation experiments were carried out
on these regions. In section 2.2, the materials and methods used to prepare the
tissue microarray map of prescribed locations, the AFM experimental setup with
its operating principle and the AFM data analysis are discussed. In section 2.3, the
tissue elasticity results are presented. This chapter is concluded in section 2.4 with
a discussion of the merits and demerits of the experimental protocol and suggested
improvements.
2.2 Material and Methods
2.2.1 Tissue Microarray (TMA) preparation and Annotation
Tissue Microarray (TMA) technology is a relatively new technique for har-
vesting tiny cylinders of tissue and arranging them on a recipient paraffin block in a
matrix-like format [16]. TMA has been widely used as a high throughput analytical
10
tool for conducting large scale molecular studies on tissue specimens. Tissue cores of
0.6 mm diameter were extracted from paraffin-embedded normal and breast cancer
tissue blocks and assembled into 4 quadrupled tissue microarray blocks using Auto
Tissue Arrayer (Beecher ATA-27). Two consecutive 4 µm slices of each TMA were
cut and fixed onto glass slides (see Fig. 2.1).
(a)
(b)
Figure 2.1: (a) Stained slide image and (b) adjacent unstained slide image c©2010
IEEE.
One of each set of consecutive slides was stained with hematoxylin and eosin
(H&E) and cover-slipped. Following the extraction on tissue cores, Virtual Mi-
croscopy (VM) technology was used to generate digital scans of the tissue segments.
VM automatically scans a specimen in a fixed high resolution and allows users to
navigate the digital specimen as if they were viewing the specimen using a tra-
ditional microscope [see Fig 2.2(a) and 2.2(b)]. Both the stained and unstained
11
tissue slides were digitized using Trestle/Zeiss MedMicro scanning system at 40x
equivalent resolution into a tiled TIFF format and uploaded onto the web server at
http://virtualscope.umdnj.edu for subsequent viewing and annotation.
The H&E slides were examined by a certified pathologist to confirm histolog-
ical validity of specimen and one pair of consecutive slides was selected based on
the fact that they contained adequate amount of breast parenchyma for the exper-
iments. The certified pathologist then annotated valid normal epithelial regions;
normal stromal regions; cancer epithelial regions and cancer stromal regions using
the online annotation tool [Fig 2.2(c)], while the labels were concealed to the person-
nel performing the AFM tests. The AFM operator was only provided information
on the regions to be indented.
Using TMA technology in the experimental design served three primary pur-
poses. First and foremost, the configuration of TMA specimen facilitated navigation
about the specimens and dictated that each sampling that was taken using the AFM
probe could be applied unambiguously to specific tissue cores by registering the
stained and unstained core images. Second, as each TMA slide can contain dozens
to hundreds of tissue cores, each potentially originating from a different patient,
TMA technology provides abundant specimen sources and experiment design possi-
bilities. Third, it has been proven that TMA technology preserves the microscopic
tissue architecture in each tissue disc for visual assessment as well as for various
qualitative and quantitative assessments [16, 20, 21, 22, 23].
12
(b) (c)
(d) (e)
Figure 2.2: An example of one tissue sample as has been used in different steps
of experiment: (a) A low resolution stained image of the tissue microarray. (b)
The VM image of an H&E stained tissue core. (c) The same core annotated by a
pathologist. (d) The unstained adjacent tissue core (VM image). Please note that
tissue distribution in this image is nearly the same as (b) but nearly impossible
to distinguish by eye. (e) Image taken from the AFM microscope during AFM
experiments c©2013 IEEE.
13
2.2.2 AFM experimental setup
Prior to AFM experiment, the unstained adjacent slide from the annotated
one was de-paraffinized with xylenes and hydrated with graded alcohols, then kept
in Phosphate buffered saline (PBS).
Top-view 
Optics
  Inverted 
microscope
Vibration 
isolation table
AFM 
head Manual XY 
positioning 
stage
Figure 2.3: AFM experimental setup c©2010 IEEE.
The Atomic Force Microscope (AFM) system consists of the AFM head and
controller (MFP-3D-BIOTM, Asylum Research) and an inverted microscope (Model:
TE2000U, Nikon, Inc) coupled to it such that the AFM head rests on the microscope
(Fig. 2.3). The whole setup is placed on a vibration isolation table (manufactured
by Herzan) and is isolated from external noise using an acoustic hood. The X-Y
14
range of the automated scanning stage is 90 µm and the Z-range is 40 µm.
The region of the slide containing the tissue cores was encircled by a hydropho-
bic barrier (Fig. 2.4) using a Pap-Pen (Invitrogen Corporation, Carlsbad, CA), and
PBS solution was intermittently added in the encircled region to ensure that the
tissue cores did not dry up. The barrier prevents the PBS solution from spilling
over to the X-Y scanning stage while performing tissue-AFM probe alignment.
Tissue hydrated in PBS
Hydrophobic
 barrier
Microscope
 Slide AFM X-Ystage
Figure 2.4: Hydrated tissue specimen on microscope slide.
After locating each array core on the slide, the annotated regions of interest
within each core were identified by closely following the annotated map [Fig. 2.2(c)]
and ∼ 5-8 indentations were carried out within the annotated region.
AFM cantilevers with pyramidal tips have been shown to produce stress con-
centration on the sample and may also potentially damage the tissue [24]. Hence,
AFM cantilevers with spherical glass tips (Novascan Technologies, spherical bead
of diameter 5 µm) were used for the AFM indentation experiments (see Fig. 2.5).
Typical probe dimensions are tabulated in Table. 2.1.
15
Chip
Cantilever
Spherical Bead
Length, 
Width, 
Thickness, 
(Radius,   )
Figure 2.5: CAD drawing of a typical AFM probe used for tissue characterization
experiments.
Table 2.1: AFM probe parameters used for characterization experiments
Typical Parameters
Length, L 130 µm
Width, W 35 µm
Thickness, t 2 µm
Bead Radius, R 5 µm
Spring Constant 4.50 N/m
2.2.3 Contact-mode operation of AFM
Fig. 2.6 shows a schematic of the contact-mode operation of the AFM. The
AFM probe is lowered vertically towards the specimen using the piezoelectric scan-
ner (also called the z-piezo), while its deflection is measured by monitoring the
deviation in laser signal reflected from the backside of the AFM probe and sensed
through a photodiode. When the probe is far away from the specimen, the vertical
z position and the deflection, d, are recorded as z1 and d1 respectively [Fig. 2.6(a)].
As the probe is lowered, it makes contact with the specimen at (zk, dk) [Fig. 2.6(b)].
The probe is further lowered till its maximum vertical displacement at (zn, dn)
16
X-Y Stage
Z-Piezo
X-Y Stage
Z-Piezo
Photodiode
X-Y Stage
Z-PiezoLaser 
Source
Photodiode
Photodiode
Mirror Mirror
Mirror
Laser 
Source
Laser 
Source
(a) (b)
(c)
3.13.23.33.43.5
−20
−15
−10
−5
0
5
10
15
20
25
30
Approach
Retract
Z-Sensor,
De
fle
ct
io
n,
(d)
Figure 2.6: Schematic of contact-mode AFM operation.
17
[Fig. 2.6(c)]. After this, the probe is retracted in the opposite vertical direction.
A 2D representation of the deflection, d, as a function of the vertical displacement,
z, constitutes an AFM force curve [Fig. 2.6(d)]. When the probe is lowered to-
wards the specimen the force curve is termed as the approach force curve. During
the probe retraction phase, the force curve is termed as retract force curve.
Given n data points in a force curve1, the net indentation in the sample is
given by the following equation [24]:
∆ = (zn − zk)− (dn − dk) (2.1)
The AFM measures the probe deflection, which needs to be related to force.
For low deflection regimes, it is possible to treat the cantilever probe as a linear
spring of spring constant kc. The value of kc is obtained by calibrating the probe’s
spring constant prior to an AFM experiment. The “thermal method” [25] was used
for calibrating the probe’s spring constant using Asylum Research software (IGOR
Pro, Wavemetrics, Inc.)2. The force, F , then can be related to the deflection by :
F = kc(dn − dk) (2.2)
Now that force F and indentation ∆ can be obtained from the AFM force
curve, all that remains is to curve-fit the force and indentation to a contact model
1Henceforth, by the term “force curve”, the approach part of the AFM force curve is implied.
In this dissertation, data analysis was not carried out on the retract portion of the force curve,
which is typically used to quantify viscoelastic properties.
2Details of the probe calibration procedure are given in chapter 4.
18
to estimate the material properties of the specimen. This is discussed in section
2.2.5.
2.2.4 Operating modes for AFM indentation
There are two modes in which the AFM force curves are acquired: deflection-
controlled mode and the indentation-controlled mode. These determine the amount
of vertical travel of the AFM probe towards (and into) the specimen.
2.2.4.1 Deflection-controlled AFM force curves
The deflection-controlled mode is considered to the be the “gold standard” for
acquiring AFM force curves, and is the most commonly used mode used for AFM
indentation experiments. In this mode, a constant deflection setpoint [(dn − d1),
see Fig. 2.7(a)] is set and the probe is lowered towards the specimen till the probe
deflection (dn) with respect to the non-contact deflection (d1) reaches the deflection
setpoint [Fig. 2.7(b)]. At this point, the probe is retracted from the specimen.
Fig. 2.7 shows the acquisition of an AFM force curve in the deflection-controlled
mode with 50 nm deflection setpoint. The deflection can be directly related to
the contact force by multiplying with the probe spring constant, consequently, the
deflection-controlled mode is also termed as force-controlled mode. Henceforth, this
mode is referred to as the force-controlled mode.
19
(a) (b)
50 nm
Figure 2.7: (a) Screenshot of MFP3D software panel used to acquire AFM force
curves in the deflection-controlled mode. (b) Acquired AFM force curve with a
deflection setpoint of 50 nm. The red and blue arrows indicate the direction of the
probe’s motion.
20
2.2.4.2 Indentation-controlled AFM force curves
For force curves acquired using the same deflection setpoint through the force-
controlled mode, specimen elasticity variations manifest as variability in indentation
(soft materials show larger indentations than stiffer materials for the same deflection
setpoint). Therefore, the force-controlled mode is not useful when a constant sample
indentation is desired.
The indentation-controlled mode is used to address this need. In order to
achieve a specific indentation in a force curve, the contact location (zk, dk) [please
see Fig. 2.6(d)] must be known. However, there is no way for the AFM to identify
the contact point in a given force curve: it can only measure the deflection, d, and
the vertical z-position, z.
To achieve a target indentation in the indentation-controlled mode, first an
AFM force curve is acquired for a very low deflection setpoint (∼ 10 nm) in the
force-controlled mode [Fig. 2.8(a,b)]. The vertical position abscissa z¯n at which
this deflection setpoint is reached is recorded [Fig. 2.8(b-inset)]. Then, a target
indentation setpoint is set [Fig. 2.8(c)] and a second force curve is acquired. In this
second force curve, the contact location zk is set to z¯n from the first force curve.
Based on this estimated contact point, the probe is lowered by the z-piezo till the
indentation [(zn− zk)− (dn−dk)] reaches the target indentation depth [Fig. 2.8(d)].
In Fig. 2.8, a deflection setpoint of 10 nm and and indentation setpoint of 250 nm
was used for the first and second force curves respectively.
While the indentation-mode appears attractive, particularly in the context of
21
Indentation = 572 - 330
               = 242 nm
330 nm
572 nm
(a)
(c) (d)
10 nm
(b)
Figure 2.8: Screenshot of cropped MFP3D software panel used to acquire AFM
force curves in the indentation-controlled mode. (a) Screenshot of 10 nm deflection-
trigger. (b) Acquired force curve with 10 nm deflection trigger. (c) Screenshot of 250
nm indentation trigger. (d) Acquired force curve with 250 nm indentation trigger.
The red and blue arrows indicate the direction of the probe’s motion.
22
analyzing force curves for mechanical characterization [shallow indentations (< 50
nm) for a 4 µm thick tissue sample) do not provide enough information about the
tissue response, while too large indentations3 necessitate complex contact models
that account for nonlinearity in the tissue response], there are significant shortcom-
ings of this approach. Firstly, the tissue specimen is indented twice in this mode,
as opposed to the force-controlled mode where only one force curve is acquired at
a given location on the tissue surface. Not only does this double the time taken
to conduct characterization experiments, but this mode also subjects the tissue to
preconditioning effects due to the effect of the first force curve. Secondly, prior
knowledge of the tissue specimen’s stiffness is necessary while acquiring the first
force curve, because even a small setpoint may result in large tissue indentation if
the tissue stiffness in the probed region is considerably less than the probe’s spring
constant. If this is the case, the deflection setpoint for the first force curve has to
lowered and the process repeated.
2.2.5 Data Analysis
The mechanical properties are obtained by processing the force-curves offline.
Two steps are needed for this: determining the contact point (zk, dk), and curve
fitting the post-contact force indentation data to a contact model.
To determine the contact point (zk, dk), a derivative-based algorithm (Fig. 2.9)
proposed in [24], was used. In this method, the slope between two deflection dat-
3The impact of large indentations on the estimated elastic modulus has been studied in detail
in Chapter 5.
23
apoints di and di+p is compared to a predefined threshold value ssetpt. If the slope
exceeds ssetpt, then the contact point is estimated to be (zkˆ, dkˆ) = (zi+p/2, di+p/2).
Based on the experimental force curves, p = 10 and ssetpt = 5− 20 nm is used.
Is
?
Yes
No
Figure 2.9: Contact estimation algorithm.
The Hertz contact model [26] was used for extracting the elastic modulus from
the acquired force curves. The contact force, F , is related to the indentation, ∆, of
an elastic half-space by:
F = 4
√
R
3
[1− ν2
E1
+ 1− ν
2
E2
]−1
∆1.5 (2.3)
where E1 and E2 are the elastic moduli of the indenter and the half-space and
ν1 and ν2 are the Poisson’s ratio of the indenter and the half-space respectively. The
contact radius, a, is given by the following relation:
a =
√
R∆ (2.4)
Hertzian theory involves the following assumptions:
24
Figure 2.10: Hertzian contact of an elastic half-space by a sphere of radius R.
1. The sphere and the half-space are homogeneous and isotropic.
2. The indentations and strains in the half-space are infinitesimal and therefore
linearized continuum theory is applicable.
3. The contact between the sphere and the half-space are considered frictionless.
4. There is no adhesion between the sphere and the indenter.
5. Both the sphere and the half-space are linearly elastic materials.
Assuming the tissue to be incompressible, ν2 can be set to 0.5. For the case of
indentation of the tissue with an probe (with an attached spherical bead), E2  E1,
hence Eqn. (2.3) can be re-written as:
F = 4E2
√
R
3(1− ν2)2
∆1.5 (2.5)
25
2.3 Results and Discussions
The AFM experiments conducted were primarily aimed at examining the dif-
ference in the mechanical properties at (a) epithelial and stromal regions of normal
breast tissue and (b) epithelial and stromal regions of cancerous breast tissue.
Three separate sets of AFM indentation studies were performed to investigate
the variations in elasticity in normal and cancerous samples. The calibrated spring
constants in each set are tabulated in Table 2.2. The experiments were conducted
in the indentation-controlled mode by setting an indentation setpoint of 250 nm.
Table 2.2: AFM probe spring constant for all experiments
Set Calibrated Spring Constant Theoretical Spring Constant
1 2.50 N/m 4.50 N/m
2 3.89 N/m 4.50 N/m
3 3.25 N/m 4.50 N/m
Figure 2.11: Histology images (20X magnification) cropped from four representative
regions probed in set 1. (a) and (b) are regions with normal epithelial (A8) and
cancerous epithelial (A13) tissue respectively, (c) and (d) are regions with normal
stromal (A22) and cancerous stromal (A12) tissue respectively c©2010 IEEE.
26
0 100 200 300 400
0
100
200
300
400
500
600
700
800
900
1000
0 100 200 300 400
0
100
200
300
400
500
600
700
800
900
1000
0 100 200 300 400 500
0
500
1000
1500
2000
2500
0 100 200 300 400 500
0
500
1000
1500
2000
2500
Indentation (nm) Indentation (nm)
Indentation (nm) Indentation (nm)
Fo
rc
e 
(n
N)
Fo
rc
e 
(n
N)
Fo
rc
e 
(n
N)
Fo
rc
e 
(n
N)
(a) (b)
(c) (d)
Hertz Fit
Experimental Data
Figure 2.12: Force-Indentation plots obtained in set 1 at (a) normal epithelial regions
(A8), (b) cancer epithelial regions (A13), (c) normal stromal regions (A22) and (d)
cancer stromal regions (A12) along with their Hertzian fits and the mean R2 values.
Inside each annotated region, roughly 5-8 force curves were taken. Note that despite
a target indentation of 250 nm, the tissue indentation varied from 150 - 500 nm.
This was primarily on account of incorrect contact estimation in the indentation-
controlled mode.
27
The force-indentation curves and the associated Hertz fits obtained by prob-
ing 5-8 locations on each of the annotations (A8, A13, A22 and A12 as shown in
Fig. 2.11) in set 1 are shown in Fig. 2.12. From the force plots, it is evident that
cancerous epithelial regions exhibit lower stiffness compared to normal epithelial
regions.
A summary of the results of all three sets of experiments is tabulated in Table
2.3. In set 1, 23 and 25 regions were probed in normal and cancerous epithelial tissue
respectively, and 12 and 18 regions were probed in normal and cancerous stromal
tissue respectively. The mean elastic modulus of the epithelial regions was computed
to be 981.80 ± 576.72 kPa and 522.79 ± 274.66 kPa for normal and cancerous tissue
respectively. The mean elastic modulus of the stromal regions was computed to be
2342.65 ± 968.38 kPa and 1335.58 ± 1087.87 kPa for normal and cancerous tissue
respectively. Using an unpaired t-test, the normal epithelial regions were found to
be significantly stiffer than cancerous epithelial regions (p = 0.006).
In set 1, the overall stiffness of 35 probed regions in normal tissue was sig-
nificantly stiffer than the overall stiffness of 53 probed regions in cancer tissue
(p=0.009). Also, the overall stiffness of 30 probed regions in stromal tissue was
significantly stiffer than the overall stiffness of 48 probed regions in epithelial tissue
(p=0.0001).
In set 2, however, only 3 cancerous stromal regions could be probed. Likewise,
in set 3, only 7 normal epithelial regions could be probed. Undersampling in the
cancerous stromal regions and normal epithelial regions in sets 2 and 3 respectively
was primarily due to (1) fewer annotated regions to sample from and (2) region
28
Table 2.3: Summary of AFM indentation experiments. Results shown as mean ± standard deviation. P-values
less than α(= 0.05) are highlighted in bold.
Epithelial Stromal
Normal Cancerous Normal Cancerous
Sampled Points 23 25 12 18
Set 1 Elastic Modulus (kPa) 981.80 ± 576.72 522.79 ± 274.66 2342.65 ± 968.38 1335.58 ± 1087.87
p-value 0.006 0.015
Sampled Points 37 13 14 3
Set 2 Elastic Modulus (kPa) 985.99 ± 1059.99 326.36 ± 128.23 2288.18 ± 1839.41 × †
p-value 0.031 × †
Sampled Points 7 21 34 26
Set 3 Elastic Modulus (kPa) × † 186.40 ± 204.21 400.03 ± 328.23 308.29 ± 241.91
p-value × † 0.236
† Elastic modulus and p-value results in these fields are excluded due to insufficient sample size arising out of undersampling in
the corresponding tissue regions.
29
localization errors arising out of manual registration methods, i.e. the AFM probings
were carried out outside the specified annotations. Hence, no conclusive inferences
could be made about the comparison between epithelial and stromal tissue response
in normal and cancerous specimens in sets 2 and 3 respectively, and therefore, the
elastic modulus values in the cancerous stromal regions in set 2 and normal epithelial
regions in set 3 have not been presented in Table 2.3.
2.4 Conclusions
In this chapter, the results of preliminary AFM characterization studies on
breast tissue specimens have been discussed. The preliminary feasibility studies in-
dicate that the AFM holds the potential to objectively distinguish between normal
and cancerous tissue. It was observed that cancer tissue exhibited significantly less
stiffness than normal tissue and that epithelial regions exhibited less stiffness than
stromal regions. It was also shown that cancerous epithelial regions were signifi-
cantly softer than normal epithelial regions. Overall results indicated that tissue
had higher stiffness as compared to previously reported measurements on cultured
cell lines or isolated single cells [9],[10]; however, this is most likely attributed to the
general tissue architecture in which the AFM studies are carried out. These results
demonstrate the feasibility of applying AFM on histologically prepared tissue sam-
ples; and this approach can potentially provide a unique pathway to gaining insight
into the biophysical changes at the onset and progression of breast cancer and other
diseases.
30
The pilot studies also highlighted certain limitations of using the current ex-
perimental AFM system used for characterizing the tissue specimens. Since the
AFM allows individual pointwise estimation of local stiffness in the tissue speci-
mens, only fraction of the tissue provided on a TMA specimen could be sampled
during each experiment setting. This is primarily due to the laborious nature of
conducting these experiments, where specimens need to be manually translated af-
ter a force curve is acquired on a tissue region. Despite the fact that the AFM user
visually registers stained and unstained tissue images to ensure AFM probings are
conducted within the specified annotated region, such manual registration methods
are not guaranteed to be completely error-free. As seen in the results in Table 2.3,
stromal and epithelial regions were significantly undersampled by the AFM in sets 2
and 3 respectively due to incorrect probing localizations. This impediment is further
aggravated when the number of annotations are limited. Reducing human interven-
tion in carrying out these tasks by automating the process of region localization
and automated indentation can increase robustness in characterizing breast tissue
specimens using AFM. This is discussed in greater detail in chapter 3.
The elastic modulus values reported in this chapter also depend on the spring
constant of the AFM probe being used. While several calibration schemes have been
proposed [25][27], all of these methods report on a certain amount of variability (∼
5−17%) [28]. Another concern with AFM indentation experiments is the sensitivity
of the estimated elastic modulus to estimation errors in the contact point [29]. These
issues are discussed in detail in chapter 4.
While the Hertz model was used to characterize the tissue specimens, it is true
31
that apriori assumptions of Herztian contact theory are not always applicable to
soft biological samples that often exhibit hyperelastic behavior [30]. However , the
focus of this preliminary work was to investigate the variations (if any) in stiffness
in normal and cancerous tissue specimens, and the Hertz model provides a rough
estimate of the stiffness of tissue specimens. Aspects of constitutive modeling of the
tissue specimens are discussed in detail chapter 5.
32
Chapter 3
Improvements to the AFM experimental protocol for tissue
characterization
3.1 Overview
In chapter 2, it was shown that using stiffness measures as a biomarker,
contact-mode AFM could be used to differentiate spatially between normal and can-
cerous histological breast tissue specimens. By visually relating stained annotated
images and unstained images acquired from AFM integrated optical microscope dur-
ing AFM sampling, manual registration of tissue regions of interest (ROI) with the
stained annotation was performed, and by using manual adjusting screws attached
to the base of the AFM stage, the tissue specimens could be aligned with the AFM
probe.
In this chapter, some of the shortcomings of using manual registration and
positioning methods are discussed in detail that were briefly alluded to at the end
of chapter 2. As opposed to pointwise tissue elasticity measures, elastic maps are
proposed as a technique to quantify the spatial distribution of tissue specimen’s
elastic modulus. The development of a semi-automated image-guided positioning
system to automate some of the positioning tasks involved in tissue characterization
is reported in this chapter. Using the developed system, it becomes possible to
33
acquire much larger experimental datasets within the same experimental setting.
3.2 Problem Statement and Proposed Solution
There were several experimental challenges that were encountered in the pre-
liminary AFM characterization studies in chapter 2. These are elaborated as follows:
• Low sample size: In the preliminary characterization studies, 78, 67 and 88
regions were sampled in experimental sets 1, 2 and 3 respectively. While the
tissue elasticity did display a trend, characterization studies for conclusive
histopathological inference would clearly require much larger datasets, poten-
tially in the order of hundreds or thousands.
• Low throughput during AFM experiments: One of the primary causes for low
sample size was the repetitive manual positioning tasks that were required
to be carried out in between two consecutive force curve acquisitions. These
included (a) moving the tissue slide so that another tissue region inside the
annotation could be sampled, (b) manually registering the stained and un-
stained images and (c) translating the tissue slide to probe annotated regions
in a different core. All these tasks carried out multiple times leads to AFM
operator fatigue, and severely restricts throughput. A flowchart showing the
various steps involved in manual tissue characterization is shown in Fig. 3.1.
• Lack of Automated Registration: As stated in the previous bullet point, man-
ual registration remains another major impediment to improving throughput
in conducting AFM based tissue characterization. The relatively large length
34
Start
Position slide underneath
AFM head such that tissue
core is in the field of view (M)
Register annotation on stained
image with unstained image
at low magnification (M)
Align annotated region
with the AFM probe (M)
Change objectives to higher
magnification and re-adjust
focusing and scene lighting (M)
Perform tissue registration
at higher magnification (M)
Change objectives to lower
magnification and re-adjust
focusing and scene lighting (M)
Is Reg-
istration
accurate
(M) ?
AFM Probing (M)
Relevant
Points in
annotation
probed
(M) ?
Manual Positioning (M)
End
NO
YES
YES
NO
Figure 3.1: Flowchart describing manual AFM indentation experiments. (M) = Manual,
(A) = Automated. 35
(b)
A1
A2
A3
(a)
(c)
A1
A2
A3
Figure 3.2: Manual registration between unstained and stained images. (a) Stained
TMA core image annotated by the pathologist as regions A1, A2 and A3. (b)
Brightfield image while performing AFM indentation on a point in annotation A3
after manually registering the stained image (c) Stained TMA core image overlaid
with the points probed during AFM experiments on the unstained slide. Observe
that annotation A2 has been incorrectly sampled by the AFM.
36
scale of the TMA cores (∼ 600 µm in diameter) implies that the whole tis-
sue core is not visible at high magnifications under the optical microscope.
Therefore it becomes necessary to zoom in on a ROI to locate annotations
desired to be phenotyped and then revert back to a coarse magnification to
get a larger perspective of the alignment of the ROI and the AFM probe (see
Fig. 3.1). Furthermore, the image of the tissue specimen acquired by the op-
tical microscope has to be registered with the stained and annotated image to
ensure the AFM probing is restricted to the annotated regions on the stained
image. Poor image contrast in the unstained slide leads to difficulty in manual
registration of the stained and unstained images and as seen in the results in
chapter 2, it might lead the AFM operator to conduct probing outside the
regions annotated by the pathologist (see Fig. 3.2).
Hypothesis: Positioning and registration automation can overcome several
of the aforementioned shortcomings of using manual AFM probe-tissue alignment
prior to tissue indentation.
The following solutions are proposed to improve the AFM experimental pro-
tocol for tissue characterization using AFM:
• Raster scanning for generating elastic maps: By acquiring a set of force curves
from points at precise spatial intervals on a tissue region of interest (ROI), a
two-dimensional map of the tissue elasticity can be constructed. Each force
curve can be processed and assigned a definite pixel of a topographic image,
where the pixel value corresponds to the elastic modulus of force curve acquired
37
at the pixel location. This resulting image is typically called an “elastic map”
[31] or a “force-volume” [18].
Elastic maps are proposed as a technique to quantify the spatial distribution
of tissue specimen’s elastic modulus. The AFM X-Y stage, which is typi-
cally employed for imaging applications, can be used for positioning the tissue
specimen at pre-specified spatial intervals.
Table 3.1: X-Y travel range of automated AFM scanners, arranged in increasing
order of their travel range.
AFM Vendor (X-Y) Stage Vendor Travel Range (µm)
Asylum Research Inc. Asylum Research, Inc. 90
Agilent 5420 Agilent 90
JPK Instruments, Inc. Tip Assisted Optics (TAO) stage 100
from Physik Instrumente, Germany
Bruker, AXS Bruker, AXS 150
• Large scanning operations with improved registration: While the commercial
AFM X-Y stages are capable of closed-loop accurate position control, they
are limited by their range of travel (∼ 90 - 150 µm) range for closed-loop
positioning (see Table 3.1). For acquiring elastic maps on small tissue ROIs,
the AFM X-Y stage might be sufficient, however, it would not be possible to
scan regions beyond a 90µm×90µm ROI. Moreover, the AFM X-Y stage still
does not resolve the registration problem.
To account for large ROI scanning operations and repeated manual registra-
tion, a semi-automated image-guided positioning system (IGPS) is proposed.
Essentially, this system uses a micromanipulator with a large X-Y travel range
38
to align the tissue specimens with the AFM probe prior to AFM indentation
using closed-loop image feedback from the camera mounted to the inverted
microscope. As described later in section 3.5, this system also addresses key
manual registration challenges.
In the subsequent sections 3.3, 3.4 and 3.5, details of the design and implementation
of the IGPS system are discussed.
3.3 AFM experimental setup with the IGPS
In the modified AFM system, an additional motorized MP-285 micromanip-
ulator (manufactured by Sutter Instruments, Novato, CA) is used for translating
the tissue slide, which is situated at the base of the microscope on the vibration
isolation table. Attached to the MP-285 is a custom-made end-effector to which the
tissue slide is mounted. The slide is placed between the AFM X-Y stage and AFM
probe. The MP-285 has a step resolution of 40 nm and a range of 2.54 cm on both
X and Y axes, and this enhanced X-Y travel range makes it possible to scan ROIs
much larger than 90µm× 90µm.
Since the forces measured during AFM indentation experiments are in the
nm-range, chatter in the end-effector or the slide can lead to incorrect mechanical
property estimation and can potentially damage the AFM probe and the tissue
sample. To eliminate vibrations in the slide, the end-effector was fabricated using
Aluminum and the slide was clamped to the end-effector using a clamping screw
(Fig. 3.3 inset). Additionally, an adhesive tape was used to firmly attach the slide
39
to the end-effector. It was observed that using this arrangement, there were no
perceptible vibrations in the slide/end-effector, as mentioned in section 3.5.
Tissue
 Slide
Adhesive 
Tape
Clamping 
Screw
AFM 
Probe
Screws for 
Specimen Alignment
Inverted
Microscope
Vibration Isolation Table
MP-285
(micromanipulator)
End-
effector
Acoustic
HoodAFM 
Head
Figure 3.3: AFM Experimental Setup with the MP-285 micromanipulator c©2013
IEEE.
3.4 Image-Guided Navigation
3.4.1 Tracking of the ROI
The presence of vision in the loop allows the AFM operator to select a certain
ROI and place the ROI underneath the AFM probe tip in an automated man-
ner. Some of the popular tracking algorithms are gradient-based methods like the
Kanade-Lucas-Tomasi (KLT) feature tracker [32] and the Scale-Invariant Feature
Transform (SIFT) [33]. However, as seen in Figs. 2.2(d) and 2.2(e) in chapter 2,
40
the unstained tissue images are characterized by low contrast and therefore unique
feature vectors are difficult to compute. In many cases, features of high contrast
(extracted from the images) have been external particles in the fluid environment,
which tend to move around randomly during the motion of the tissue slide under-
neath the AFM probe (see Fig. 3.4).
(a) (b)
Figure 3.4: Extracted features with strong intensity gradients. Observe that promi-
nent features are impurities floating in the liquid environment of the sample.
As a result, normalized cross-correlation based template matching algorithm
[34] is used to ensure robust tracking of the ROI during the motion of the slide (see
Fig. 3.5). Though the applicability of such methods are conditional upon uniform
scene lighting and in-plane translation without rotation, these conditions could be
enforced by: (1) manually adjusting the external lighting and (2) ensuring that there
is no relative motion in the various interconnected parts of the end-effector.
The implementation of the tracking algorithm is given as follows:
For an image, It = {I(px, py, t)|0 ≤ px ≤ R, 0 ≤ py ≤ C}, and a ROI,
T = {T (px, py, t)|0 ≤ px ≤ r, 0 ≤ py ≤ c}, where It and T are the image and the
ROI respectively at the tth frame, the estimated position of the ROI at the t+ 1th
41
(a) (b)
Figure 3.5: Normalized cross correlation algorithm used for tracking at two time
instants (a) and (b). The top images correspond to It = I(px, py, t) with the ROI
marked with a white box, while the bottom ones are the result images R, where
with the estimated ROI position (pˆx,t+1, pˆy,t+1) appear as spots of high intensity.
frame is given by [34]:
(pˆx,t+1, pˆy,t+1) = argmaxpx,py R (3.1)
where R, the result image is given by:
R =
∑
p′x,p′y
[T (p′x, p
′
y)It+1(px + p
′
x, py + p
′
y)]
√∑
p′x,p′y
T (p′x, p
′
y)2
∑
p′x,p′y
It+1(px + p
′
x, py + p
′
y)2
(3.2)
The numerator of the expression in Eqn.(3.2) indicates the correlation of the image
and the ROI, while the denominator is the normalizing term to ensure that general
lighting differences in both the image and the ROI do not affect the tracking algo-
rithm significantly. The tracking algorithm was implemented in Visual C++ using
42
OpenCV libraries [35].
3.4.2 Control Scheme
The control law in discrete state space form is given by:
X(k + 1) = X(k) + u(k) (3.3)
Y (k) = SRX(k) (3.4)
where
S =


sx 0
0 sy

 ; R =


cos(θ0) −sin(θ0)
sin(θ0) cos(θ0)

 (3.5)
Y =
[
px py
]T
and X =
[
x y
]T
represent the image and manipulator
frame coordinates respectively. R is the rotation matrix between the manipulator
and the image frame while S is a scaling matrix relating the two frames. The
parameters of R and S are estimated in a pre-experiment calibration step using a
cover-slip with uniform grids marked at 100 µm separation, as shown in Fig. 3.6.
u(k) is the control input in the manipulator frame.
To estimate the control input, a gradient-descent based approach [36] is used,
where the control input u(k) is based on the gradient of the function:
F [X(k)] = [X(k)−X tip]T [X(k)−X tip] (3.6)
43
Figure 3.6: Calibration of the scaling matrix S using a cover-slip with uniform grids
at 100 µm.
The control input u(k) is therefore given as:
u(k) = −γ(k)∇F [X(k)] = −2γ(k)[X(k)−X tip] (3.7)
where the step size γ(k) is given by:
γ(k) =



γ0
2||X(k)−X tip||
, if ||X(k)−X tip|| > 
0, if ||X(k)−X tip|| ≤ . (3.8)
γ0 is a constant that determines the magnitude of the incremental travel of
the ROI towards the probe tip X tip. The positioning is terminated once the ROI is
within a radius of  from the probe tip, where  is a preset parameter.
It is vital to ensure that parts of the tracked ROI are not occluded by the
AFM probe during positioning since the template matching algorithm requires the
entire ROI to be visible in the image space. As a result, the estimated position of
44
the slide in image coordinates, Y , is selected to lie on the upper-corner of leading
edge of the ROI, as shown in Fig. 3.7(b).
3.5 ROI positioning protocol
Tracking
ROI
Probing
ROI
Coarse
ROI
(b)
(c) (d)
(a)
Figure 3.7: Tracking protocol for probe-ROI alignment at 10X and 20X. (a) Stained
image of tissue core annotated by the pathologist at low resolution (10X). (b) Bright-
field image from the AFM-optical microscope. A coarse ROI is selected to match
one of the annotated regions in (a). (c) Positioning at 10X. (d) Registration of the
coarse ROI at 20X c©2013 IEEE.
Using the tracking algorithm and the control law described in section 3.4, the
following protocol is implemented for alignment of the tissue ROI and the AFM
probe tip (see Fig. 3.7).
45
• A coarse ROI is selected at low magnification (m1 = 10x), after visually cor-
relating the annotated regions from the stained image [Fig. 3.7(a)] and the
brightfield image from the AFM-optical microscope [Fig. 3.7(b)]. The AFM
probe tip is also selected visually at the same magnification [Fig. 3.7(b)].
• Based on the tracking algorithm and control law discussed, the ROI is posi-
tioned within an error of , chosen to be 2 µm [Fig. 3.7(c)].
• The objectives, lighting and focusing are altered manually at higher magnifi-
cation (m2 = 20x) to ensure that the probe tip and the ROI is in focus and
the scene is uniformly lit. The AFM tip is then selected visually by the user
at high magnification, m2. Since the probe tip is stationary, it serves as a ref-
erence for the coarse ROI, which is recreated at a distance m2/m1 from the
tip at m2 magnification [Fig. 3.7(d)]. This allows registration of the same ROI
across multiple magnifications. At m2 magnification, finer details are visible
to the user and a part of the recreated ROI is selected, called the probing ROI,
which is probed by the AFM. The remaining part of the recreated ROI serves
as the tracking ROI. The probing ROI is then sampled in a raster fashion,
while the tracking ROI is used to provide image feedback [Fig. 3.7(d)].
In the IGPS, the user has to visually register the stained and unstained images
only twice during the course of the experiment: first, during the selection of the
Coarse ROI at low magnification and second, during the selection of the probing
ROI at higher magnifications (see Fig. 3.8). This is a marked improvement over
manual registration procedures, where the AFM user has to visually register the
46
Start
Calibrate the elements
of S and R matrices (M)
Position slide underneath
AFM head such that tissue
core is in the field of view (M)
Locate the AFM probe tip
and register coarse ROI with
annotation on stained image
at low magnification (M)
Position coarse ROI so that
the leading edge is within an
error  from the AFM tip (A)
Change objectives to higher
magnification and re-adjust
focusing and scene lighting (M)
Perform tissue registration
at higher magnification (A)
Select the probing ROI (M)
Raster Positioning (A)
AFM Probing (M)
Relevant
Points in
probing
ROI
probed
(A) ?
End
YES
NO
Figure 3.8: Flowchart describing AFM indentation experiments using the image-
guided positioning system (IGPS). (M) = Manual, (A) = Automated.
47
annotated ROI on the unstained image periodically to ensure that the AFM probings
are carried out within the specified annotations (see Fig. 3.1).
A representative AFM force curve obtained from samples mounted on the
IGPS is shown in Fig. 3.9(a) and Fig. 3.9(b). Minimal fluctuations in the deflection
data is indicative of negligible vibrations introduced due to the end-effector design.
3.13.23.33.43.5
−20
−10
0
10
20
30
Approach
Retract
0 50 100 150 200
0
20
40
60
80
100
120
Fo
rc
e
Indentation
De
fle
ct
io
n
Z-Sensor
Experiment
Fit
(a) (b)
0 20 40 60 80
−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Error in X
Error in Y
Step Positions
Er
ro
r
(c)
Figure 3.9: (a) Representative AFM force curve on a tissue sample mounted on the
end-effector. (b) Corresponding force-indentation curve overlaid with the Hertzian
fit and (c) Tracking Performance c©2013 IEEE.
48
The positioning accuracy is demonstrated in Fig. 3.9(c), where the slide is
translated from its original position (262,226) to Y tip = (163, 362) in the image
space. The positioning errors at the completion of alignment at 10x and 20x were
1.6 µm and 0.8 µm respectively, which are acceptable alignments errors.
3.6 Results and Discussion
A combination of the AFM X-Y stage and the IGPS was used to obtain elastic
maps from the tissue specimens. In the case of the AFM X-Y stage, the annotations
on the stained slide were manually registered with the unstained slide and a 60µm×
60µm ROI inside the annotation was raster-scanned. In the case of the IGPS, the
protocol given in section 3.5 was followed to define the probing ROI on which raster-
scanning was performed.
Due to the relative ease of the force-controlled mode, this mode was used to
acquire the force curves and a deflection setpoint of 50 nm was used. Using the
contact point estimation scheme [24] and the Hertz contact model [26] described in
chapter 2 to fit the experimental data, the AFM force curves were post-processed
to obtain the elastic maps.
3.6.1 Characterization of Epithelium-Stroma boundary
Due to the sharp differences in elasticity in the epithelial and stromal tissue
regions in both normal and cancer specimens, the boundaries of these two tissue
types were identified and raster-scanned.
49
kPa
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
kPa
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
kPa
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
(a) (b)
(c) (d)
(e) (f)
Figure 3.10: Elastic maps on normal tissue cores along with their correponding
brightfield images.
50
kPa
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
kPa
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
kPa
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
(a) (b)
(c) (d)
(f)(e)
Figure 3.11: Elastic maps on cancer tissue cores along with their correponding
brightfield images.
51
Raster-scanning results on the epithelium-stromal tissue boundaries on normal
and cancer specimens are shown in Figs. 3.10 and 3.11 respectively. The first
column in Figs. 3.10 and 3.11 correspond to the brightfield image of the probed
region on the tissue core, while the second column represents the corresponding
elastic map. The first row images in Figs. 3.10 and 3.11 correspond to elastic maps
on probing ROIs (115µm× 160µm) and (90µm× 130µm) on a normal and cancer
core respectively, acquired by the IGPS and sampled at 5µm spatial intervals. The
second and third row images in Fig. 3.10 and Fig. 3.11 correspond to AFM X-Y
stage acquired 60µm× 60µm elastic maps sampled at 2.5µm spatial intervals.
In all the elastic maps, a clear delineation between the two tissue regions is
clearly observed: the softer regions correspond to epithelial tissue, while the stiffer
regions are the stroma. These results validate the observations in chapter 2, where
epithelial tissue in both normal and cancerous specimens were more deformable
compared to stromal tissue.
7.4 7.6 7.8 8 8.2 8.4
160
180
200
220
240
260
280
300
320
340
Z-Sensor
De
fle
ct
io
n 
Figure 3.12: Uninterpretatble force curve with no prominent non-contact regime.
It should be noted that during the raster-scanning, some force curves were
52
acquired where no defined contact point existed, as shown in Fig. 3.12. These
types of force curves are typically acquired when the tissue specimen has started
translating before the probe has completely retracted from the specimen surface
[37]. These uninterpretable force curves were not considered in constructing the
elastic map, and pixels corresponding to these force curves were replaced by blank
white spaces, as seen in Figs. 3.10(d), Fig. 3.10(f) and Fig. 3.11(d).
3.6.2 Characterization of patient specific tissue specimens
The preliminary characterization experiments in chapter 2 were conducted on
normal and cancer tissue cores compiled from multiple patients. In this section,
patient-specific tissue characterization results are reported, primarily to investigate
if there was any tendency of the tissue response depending on the patient from
whom the cores were extracted.
Table 3.2: AFM experimental setting parameters.
Set Patients Calibrated Spring Constant Theoretical Spring Constant
1 A,B,C 0.596 N/m 0.95 N/m
2 D,F,G,H 0.589 N/m 0.95 N/m
Two sets of experiment were conducted to investigate the tissue response from
7 patients (details in Table 3.2). Patient-specific results are tabulated in Table
3.3. On observation, it can be seen that that stromal tissue is stiffer compared to
epithelial tissue, as observed previously during raster scans at the epithelium-stroma
interface.
53
Table 3.3: Summary of AFM raster-scanning results on patient specific tissue specimens. Results expressed as
mean ± standard deviation.
Epithelial Stromal
Normal Cancerous Normal Cancerous
Patient A Sampled Points × 2285 × 537Elastic Modulus (kPa) × 419.65 ± 312.41 × 1149.10 ± 675.10
Patient B Sampled Points × 1144 × ×Elastic Modulus (kPa) × 461.86 ± 264.09 × ×
Patient C Sampled Points 576 × 565 ×Elastic Modulus (kPa) 512.34 ± 251.53 × 902.73 ± 524.86 ×
Patient D Sampled Points × 1152 × ×Elastic Modulus (kPa) × 359.49 ± 236.16 × ×
Patient E Sampled Points × 2304 × ×Elastic Modulus (kPa) × 152.20 ± 111.46 × ×
Patient F Sampled Points × × 576 ×Elastic Modulus (kPa) × × 714.82 ± 348.07 ×
Patient G Sampled Points 3452 × × ×Elastic Modulus (kPa) 240.61 ± 144.15 × × ×
54
It was also observed that cancer epithelial and normal epithelial tissue response
was not uniform across patients. While patients A and B show similar mean elastic
modulus in cancer epithelial regions, patients D and E show vastly dissimilar elastic
modulus values. Cancer epithelial regions in Patient E showed significantly lower
tissue response compared to normal epithelial regions in Patient G, however, when
the cancer epithelial response from Patient D was included (conducted in the same
experimental setting), the hypothesis of reduction in elastic modulus with cancer
progression could no longer be considered valid.
These results clearly indicate that normal and cancer tissue extracted from the
same patients are necessary to draw meaningful conclusions about tissue property
alterations with cancer progression. Due to time constraints and unavailability of
paired tissue specimens from the same patient, tissue response from the same patient
could not be investigated, and this remains a subject of future work.
3.7 Conclusions
In this chapter, several improvements to the original AFM experimental pro-
tocol was proposed. The major contributions in this chapter are as follows:
1. Improved Sample Size: Using a combination of the AFM X-Y stage and an
image-guided positioning system (IGPS), sample size (upto two orders in mag-
nitude) could be significantly increased, thus confirming the research hypoth-
esis of this chapter. The AFM X-Y stage is employed to perform small scans,
while the IGPS allows scans on much larger ROIs. This is especially relevant to
55
the study of relatively large histopathological breast tissue microarrays, where
malignancies of interest are not easily isolated within a small tissue region.
2. Improved Registration: The use of the image-guided positioning system allows
significant reduction in manual registration of the stained and unstained im-
ages. The user requires to perform manual registration only twice during the
course of the AFM experiments, as opposed to reviewing the position of the
ROI on the unstained image with respect to the annotated image periodically
in manual positioning operations.
There are, however, certain limitations with the IGPS, which are elaborated
as follows:
1. User interventions in microscope operations required: Though manual posi-
tioning operations could be significantly reduced, the AFM user is still re-
quired to shift objectives, alter the objective focal length and adjust the scene
lighting manually. It is possible to automate these tasks, however, these au-
tomated options are not currently available in the inverted microscope (Nikon
TE2000U).
2. Absence of automated registration between the stained and unstained images:
Further user intervention could be reduced by performing automated registra-
tion between the annotated regions in the stained and the unstained tissue
images. This would require performing multimodal registration since the two
images are acquired by different scanning systems in dissimilar lighting envi-
ronments.
56
Larger datasets implied that stronger conclusions could be drawn from the
AFM measurements. While delineations at the boundaries of epithelium and stroma
tissue could be clearly seen, it was also observed that the tissue properties also
depend on the patient from whom the tissue cores were extracted. It is therefore
proposed that normal and cancerous tissue pairs from the same patient be extracted
and characterized by the AFM.
57
Chapter 4
Probabilistic Estimation of the Elastic Modulus
4.1 Overview
In chapter 3, it was demonstrated that by using an image-guided positioning
system with a range much larger than commercial AFM stages, considerable hu-
man intervention could be reduced in AFM tissue characterization experiments. In
addition to reducing the time taken to conduct the experiments, large-scale AFM
indentation data acquisition was shown possible. To exploit the full potential of the
image-guided automated positioning system, it is necessary to develop robust com-
putational methods to automate the mechanical property extraction from large-scale
raw AFM force curves.
As discussed in chapter 2, an accurate material characterization procedure
from the AFM experimental data entails unequivocal determination of the following
experimental parameters: the contact point on the AFM force curve and the spring
constant of the AFM probe. A derivative-based algorithm was used to determine
the contact point [24], and the thermal method was used to calibrate the AFM
probe spring constant [25] to obtain the results in chapter 2. The challenges in
accurate determination of the contact point and the probe spring constant were
briefly alluded to at the conclusion of chapter 2.
In this chapter, the sensitivity of the estimated mechanical property to uncer-
58
tainties in these two experimental parameters are discussed in detail. A probabilistic
mathematical model based on Bayesian analysis is proposed to estimate the mechan-
ical property after accounting for observed uncertainties in the contact point and
the probe spring constant. In section 4.2, causes for uncertainty in the experimental
parameters are reviewed. In section 4.3, the details of the mathematical model are
described. In section 4.4, results of implementation of the model on simulated and
experimental AFM force curves are shown. This chapter closes with discussions of
the proposed model with its merits and demerits in section 4.5.
4.2 Problem Statement and Proposed Solution
4.2.1 Uncertainty in the contact Point
The determination of the contact point is a vital component of the mechanical
property extraction process because it determines the amount of indentation caused
due to the vertical travel of the AFM probe. A representative AFM force curve
on a fixed mouse Embryonic Stem Cell (mESC), is shown in Fig. 4.1. It is evident
from Fig. 4.1 that the force curve displays a smooth transition from the non-contact
regime to the contact regime. This is typical of compliant biological samples like
cells and tissue, and this smooth nature of the transition leads to considerably
difficulty in ascertaining the contact point. On stiff substrates in air, the AFM force
curve while approaching the sample shows a prominent attractive region (primarily
contributed by the van der Waal’s forces and capillary forces [38]), which causes
the first derivative of the force curve to change sign when the contact forces take
59
over [39]. However, this approach is invalidated for compliant specimens like cells
and tissue in liquid environments. As seen in the Fig. 4.1, there is no perceptible
attractive region in the approach curve while indenting biological specimens.
11.5 12 12.5 13 13.5 14 14.5
−300
−250
−200
−150
−100
−50
0
(a)
De
fle
ct
io
n 
12.8 12.9 13 13.1 13.2 13.3 13.4
−300
−295
−290
−285
−280
−275
−270
−265
−260
(b)
De
fle
ct
io
n 
Z-position 
1 2
3
4
5
6
7
Figure 4.1: Representative AFM force curve on mESC, with candidate contact
points shown in red blocks: (a) Whole AFM force curve and (b) transition region
from non-contact to contact regime.
The effect of uncertainty in the contact point on the elastic modulus computed
60
using a Hertzian fit [24] is shown in Table 4.1. Over a wide range of candidate contact
points, the elastic modulus varies by a factor of 42% [expressed as (min - max)/ave
].
Table 4.1: Variation of elastic modulus with change in contact point.
Contact Point 1 2 3 4 5 6 7
Elastic Modulus (kPa) 13.28 14.29 15.35 16.50 17.72 19.03 20.52
Over the years, researchers have proposed various strategies to estimate the
point of contact in AFM force curves. Most contact estimation algorithms can be
broadly divided into three categories : purely fit-based, purely derivative based and
a combination of fit and derivative based algorithms. Following is a brief description
of the existing algorithms in literature, together with their strengths and weaknesses.
• Purely Derivative-based Approaches
– Pillarisetti et al. [24] and Nyland and Maughan [40] assumed contact
when the first derivative of the deflection with respect to the Z-position
exceeded a threshold.
– Radermacher [41] detected the contact point as the point of discontinuity
in the first derivative of the raw force curve.
– Advantages: These methods are independent of the contact-model used
to extract the material properties.
– Disadvantages: Computing derivatives in the presence of noisy data could
lead to incorrect estimates. Moreover, smoothness of the transition re-
61
gion might preclude any prominent changes in the first derivative of the
deflection. In addition, selection of the threshold parameter to detect the
change in the deflection derivative [24],[40] is largely subjective and is
prone to noise in the non-contact regime.
• Derivative and fit-based Approaches
– Jaasma et al. [39] used a piecewise quadratic function to fit the deflection
data in the force curve, and estimated the contact point as the point
where the extrapolated first derivative of the fitting function exhibited
zero-crossings.
– Charras et al. [42] and Radermacher et al. [43] solved for the contact
point and the elastic modulus simultaneously by taking two points on the
force curve, whose existence in the contact regime was beyond doubt.
– Dimitriadis et al. [44] performed a unconstrained sequential search for
the contact point by fitting the force curve to a modified Hertz contact
model, and selected the point that led to the best fit as the contact point.
– Lin et al. [29] used high-order derivatives to obtain a truncated dataset
where linear elasticity theory was valid. Following this, a sequential
search for the contact point was performed by fitting all the points to
the right of the candidate contact point. The datapoint that led to the
best fit was chosen as the contact point. Crick and Yin [45] used a secant
method to truncate the AFM force datasets to ensure applicability of
linear elasticity theory, following which a fitting protocol similar to [29]
62
was used to estimate the contact location.
– Advantages: These methods, when used on smoothed datasets, are inde-
pendent of noise in the deflection data.
– Disadvantages: Fitting (truncated) datasets to pre-established contact
models makes implicit assumptions on the nature of material being in-
dented. In many cases, the contact model is not beforehand.
• Purely fit-based Approaches
– Costa et al. [1] used a bi-domain polynomial (BDP) fit (linear pre-contact
and Hertzian contact regime) for the whole dataset. The estimated con-
tact point was the point that resulted in the best fit for both regimes.
– Polyakov et al. [46] performed segmentation of the AFM force curve
into electrostatic and contact regions using a least-squares optimization
approach.
– Advantages: These methods do not require truncation of the datasets,
which is often a subjective exercise. Also, absence of derivative compu-
tations eliminates the need of smoothing the raw datasets.
– Disadvantages: Fitting datasets to pre-established contact models makes
implicit assumptions on the nature of material being indented. In many
cases, the contact model is not beforehand.
Research has also been directed towards developing methods that obviate the
need to accurately determine the contact point. Al-Hassan et al. [47] developed an
63
approach termed as Force Integration to Equal Limits (FIEL) to obtain qualitative
maps of relative elasticity in a biological sample. Briefly, the FIEL approach involves
integrating the raw AFM force curves to the same limits (i.e. maximum deflection).
The authors reported that varying the contact point by a few datapoints resulted
in very small changes in the computed integral. The stiffness at different spatial
locations was then expressed an a non-dimensional parameter, normalized by the
integral at a given spatial location.
A relatively new approach to estimate the contact point was recently proposed
by Rudoy et al. [48], which involves a Bayesian Changepoint model to obtain pos-
terior distributions of the contact point along with the elastic modulus. As opposed
to the other approaches stated before which are primarily point-based, a Bayesian
formulation generates distributions that makes it possible to investigate the vari-
ability in the elastic modulus caused due to the contact uncertainty. The Bayesian
formulation proposed in [48] can be looked upon as a probabilistic extension to the
bi-domain polynomial (BDP) fitting approach proposed by Costa et al. in [1], which
is a maximum-likelihood based approach to estimate the contact point.
4.2.2 Uncertainty in AFM Probe Calibration
As stated previously in chapter 2, the deflection data, d, is related to the
elastic modulus using the Hertz model as shown below:
F = kc(dn − dk) =
4E
√
R
3(1− ν2)∆
1.5 (4.1)
64
where k is the contact point index in a dataset of n datapoints, kc is the AFM
spring constant, E is the elastic modulus and ν is the Poisson’s ratio of the material
being indented. Since the elastic modulus E scales with the spring constant kc, any
uncertainty in kc leads to a proportional change in E.
Due to potential inaccuracies during microfabrication, AFM probes are typi-
cally calibrated prior to their use. Various calibration methods have been proposed
by researchers. Some of the most prominent ones are:
1. Sader’s Method [27]: This method uses the probe dimensions, fundamental
resonance frequency and the quality factor of the noise spectrum of a cantilever
vibrating due to the thermal energy of the surroundings to calibrate the spring
constant of the probe.
2. Thermal Method [25]: Similar to Sader’s Method, the thermal method utilizes
equipartition theorem to estimate the probe spring constant of a cantilever
vibrating due to the thermal energy of the surroundings.
3. Cleveland Method [49]: The Cleveland method determines the probe spring
constant from the decrease in resonant frequency resulting from the attach-
ment of known masses to the probe.
4. Reference cantilever based methods [50]: This method involves deforming a
probe of unknown spring constant against a known (reference) probe, and then
measuring the deformations to estimate the spring constant of the unknown
probe.
65
Of the aforementioned calibration protocols, the first three are most com-
monly used to estimate the spring constant. Of these, the thermal method [25] has
emerged as a preferred choice of many researchers over (a) the Sader method [27]
because it requires no prior knowledge of accurate dimensions of the probe and (b)
the Cleveland method [49] because of its relative ease of use, as opposed to the
cumbersome procedure of mounting masses on AFM probes [49]. Consequently, the
thermal method is used to calibrate the AFM probe spring constant using Asylum
Research software (IGOR Pro, Wavemetrics, Inc.)
The thermal method involves determination of (1) the sensitivity of the pho-
todiode, commercially called Inverse Optical Lever Sensitivity (InvOLS) [nm/v] and
(2) fundamental resonant frequency of the probe. Following is a step-by-step process
to obtain the spring constant:
• An AFM force curve is obtained on a tissue-free hard glass surface of the
microscope slide [Fig 4.2(a)].
• The slope of the approach force curve is computed.
• The InvOLS is computed from the slope of the force curve, after incorporating
correction factors that accounts for the approximation of AFM probes as ideal
springs and finite laser spot size on the probe end [51].
• The probe is withdrawn from the surface of the microscope slide and is made
to vibrate in response to the thermal noise in the fluid surroundings.
• A Lorentzian function is then fitted to the thermal spectrum after visually
66
9.29.49.69.810
−20
0
20
40
60
80
100
120
0 50 100 150
1.5
2
2.5
3
3.5
4
x 10−13
Approach
Retract
Slope (related to 
InvOLS)
Exp. 
Fit
Z-position 
De
lfe
ct
io
n 
(n
m
)
Am
pl
itu
de
 
Fundamental
Freq. (kHz)
(a)
(b)
Figure 4.2: AFM Probe Calibration using the thermal method. (a) An AFM force
curve on a tissue-free hard glass surface and (b) thermal spectrum of the probe in
PBS solution. The rectangular probe of Table 4.2 (rated compliance of 0.222 m/N)
was used to generate these results c©2014 IEEE.
67
locating the position of the fundamental frequency [Fig 4.2(b)]. The integral
of the Lorentzian, P , is then computed.
• The spring constant is finally computed as:
kc =
kBT
P (InvOLS)2 (4.2)
While the calibration protocol discussed is fairly straightforward, interpreting
the calibration results is a challenging exercise. Most calibration schemes report
variability between 5-17%[52]. Carrying out calibration experiments in liquid envi-
ronments add to further complexity in accurate determination of the probe’s spring
constant. Calibrated spring constants have shown large variations depending on the
viscosity of the liquid medium. Pirzer and Hugel [53] reported an error of 25% in 4M
phosphate solutions and upto 100% in highly viscous solutions. It was observed that
improving the thermal spectrum fitting function to estimate the spring constants
led to reduced errors in the calibrated stiffnesses.
A relatively new technique of AFM probe calibration is the use of laser doppler
vibrometry (LDV) based interferometric methods [54] to measure the actual vertical
displacement of the probe tip, as opposed to the optical detection methods used in
standard AFMs that measure the angular deflection of the tip end. LDV-integrated
AFMs have recently gained a lot of popularity for dynamic AFM methods in liquids
[55]. Gates el at [28] recently reported calibration errors within 2%, compared to
around 5% at best for the thermal method of calibrating AFM probes; however,
68
Table 4.2: Summary of AFM probe calibration results using the thermal method.
Results expressed as mean ± standard deviation.
Probe Rated Comp- InvOLS Fundamental Calibrated
Geometry liance (m/N) (nm/V) Freq. (kHz) Compliance (m/N)
V-Shaped 16.67 87.72 ± 3.66 3.77 ± 0.13 16.38 ± 1.63
(Probe 1)
Rectangular 0.222 178.36 ± 1.88 62.68 ± 2.75 0.328 ± 0.032
(Probe 2)
LDV-integrated AFMs have yet to be commercialized, and the thermal method is
used for calibrating the AFM probes in this dissertation.
It is worth mentioning here that henceforth, the probe spring constant
is referred to by its inverse, i.e. probe compliance, s = 1/kc. Such a
transformation is a necessary step for expressing the uncertainty in the AFM probe
calibrations as an error-in-variables (EIV) regression problem, as seen later in this
chapter.
A summary of 20 calibration experiments using the thermal method for two
probes in PBS solution prior to and after the completion of the AFM experiment are
shown in Table 4.2. As seen from the thermal noise spectrum in [Fig 4.2(b)], the
location of the fundamental frequency is not sharp due to the low Q-factor of the
surrounding fluid. As a result, proper fitting of the noise spectrum to a Lorentzian
is challenging. From Table 4.2, it can be observed that the probe compliance shows
10% errors, which are in good agreement with previously reported variability with
the thermal method [52].
The aforementioned factors of uncertainty clearly pose a serious hindrance in
obtaining accurate estimates of the mechanical properties of specimens subject to
69
AFM indentation. This need is further underscored in light of fact that mechanical
characterization results on biomaterials such as breast tissue specimens are quite
often followed by statistical hypothesis testing for inference purposes, for example,
t-tests [10]. Pointwise estimation of the mechanical properties without accounting
for the underlying uncertainty can lead to erroneous conclusions, especially when the
computed p-values for the chosen hypothesis test are close to the level of significance,
α. Indeed, hypothesis testing methods that incorporate interval uncertainty [56]
reduce the possibility of inference errors compared to conventional hypothesis testing
techniques which assume that the exact values of the estimates are known. An
accurate inference test, however, has to be preceded by an approach that quantifies
the elastic modulus with its associated interval uncertainty in a robust manner [57].
Error quantification in AFM studies has been recently shown by [58], where
studies were conducted in air to measure the elastic modulus of relatively stiff sub-
strates like cellulose nanocrystals. However, the effect of contact uncertainty was
not investigated, which is considerably lesser in air compared to liquid.
Hypothesis: Accuracy of the estimated elastic modulus can be significantly
improved by developing robust uncertainty quantification techniques using Bayesian
methods.
In this chapter, a rigorous statistical approach using an integrated error-in-
variables (EIV) based Bayesian Changepoint formulation is proposed to estimate
the elastic modulus of biological samples as a function of both the contact point
and the AFM probe compliance variability. The probe compliance variability is
modeled using the EIV formulation, which is typically used in cases where the in-
70
dependent variables in a regression model are observed with errors. Use of Markov
Chain Monte Carlo (MCMC) methods makes it possible to factorize the joint pos-
terior distribution and obtain samples from the marginal posterior distributions.
The proposed algorithm is validated on simulated data and implemented on AFM
datasets obtained from indentation studies on breast tissue specimens. Finally, a
sensitivity analysis is carried out to monitor the performance of the algorithm to a
wide range of probe compliance values.
4.3 Bayesian Analysis
4.3.1 Transformation of raw data
The AFM raw data comprises the probe deflection (d) and z-position (z),
which are processed offline to estimate the properties of the specimen being studied.
Using the same notation used in section 2.2.3 of chapter 2, the contact point is given
as (zk, dk), where k ∈ (1, n) is the unknown index in the dataset which indicates the
transition from non-contact to the contact regime.
The following transformation is used:
δi = zi − di (4.3)
Physically, δ represents the tip-sample distance prior to contact and the spec-
imen indentaion post-contact. Using Eqn. (2.1), the indentation in the sample can
71
3 3.5 4 4.5
−20
−10
0
10
20
30
40
50
De
fle
ct
io
n
Position
Z-position
Figure 4.3: Representative AFM force curve of d v/s z and d v/s δ. Please note that
the d v/s δ curve is displaced 10nm on the deflection axis for ease of visualization
c©2014 IEEE.
be written as:
∆ = δn − δk (4.4)
The experimental data can now be related to an explicit force-indention rela-
tionship by the following expression:
kc(di − dk) = f(E, ν, δi − δk) (4.5)
where kc is the probe spring constant, k + 1 ≤ i ≤ n, f(E, ν, δi − δk) is the
contact model, E is the elastic modulus and ν is the Poisson’s ratio.
Two contact models are used in the EIV-Bayesian Changepoint analysis and
are given as follows:
1. Hertz contact model [26], as discussed in chapter 2.
72
2. Long’s contact model [59], for finite indentations on a thin hyperelastic speci-
men following a neo-Hookean constitutive material model.
Assuming incompressibility in the biological specimens, ν can be set to 0.5.
Both Hertz and Long’s contact models are linear in the elastic modulus, E, and can
be written as:
kc(di − dk) =
16E
√
R
9 fi (4.6)
where
fi =



(δi − δk)1.5 for f ≡ Hertz
[
(δi − δk)1.5 + 1.15
√γ(δi − δk)2 + α1γ
√γ(δi − δk)3 for f ≡ Long’s Model
+α2γ3(δi − δk)4.5
]/[
1 + 2.3γ√γ(δi − δk)1.5
]
(4.7)
The constants α1, α2 and γ are given by:
α1 = 10.05− 0.63
√
h/R(3.10 + h2/R2) (4.8a)
α2 = 4.80− 4.23h2/R2 (4.8b)
γ = R/h2 (4.8c)
where k + 1 ≤ i ≤ n. R is the of radius of the spherical bead (= 2.5 µm)
and h is the thickness of the tissue sample. Long’s force-indentation relationship is
valid in the regime ∆/h ≤ min(0.6, R/h) and 0.3 ≤ R/h ≤ 12.7. The expressions in
Eqns. (4.8a) and (4.8b) assume frictionless contact between the AFM probe and the
73
specimen, which is a reasonable assumption since the tissue specimens are hydrated
by PBS solution during the AFM experiments.
While using the Hertz model to study cells of radius Rcell, R refers to the
effective radius, given by:
R =
√
RbeadRcell
Rbead +Rcell
(4.9)
An alternate form of Eqn. (4.6) is obtained by using the compliance of the
probe, s = 1/kc, in the right-hand-side of the Eqn. (4.6). Therefore, one can write:
(di − dk) =
16E
√
R
9 sfi (4.10)
Such a transformation allows us to relate the deflection, d, to the transformed
variable δ, and to incorporate the probe compliance, s, into the EIV model, as shown
later in this section.
In the non-contact regime, the deflection, d, can be modeled as a linear function
of z [1] , i.e. d = az + b, where a and b are arbitrary constants. However, using a
simple rearrangement of the coefficients a and b, it can be readily shown that d also
varies linearly with the transformed variable δ in the non-contact regime:
d = a1− aδ +
b
1− a (4.11)
The two regime regression model relating the deflection d to δ can therefore
74
be written as:
di =



[
1 δi
]
β1 + 1 if i ≤ k
[
1 sfi
]
β2 + 2 if k + 1 ≤ i ≤ n
(4.12)
where j ∼ N(0, σj2); j = 1, 2 are independent and identically distributed
(i.i.d) normal random variables and fi is given by Eqn. (4.7). In general, σ21 6= σ22 ,
since σ21 results from the viscous interactions between the probe and the liquid
medium, while σ22 depends primarily on the probe-sample frictional forces [48]. βj =
[
βj1 βj2
]T
; j = 1, 2 are the regression coefficients in the non-contact and contact
regime respectively, and β22 encapsulates the elastic modulus E.
4.3.2 Classical Error-In-Variables Model
In a standard linear regression model with a single independent covariate ξ,
the data comprises n observations of the response variable y and ξ. These are related
by the data error model, and is given by:
yi = β0 + β1ξi + i ; i = 1, 2, 3, . . . , n. (4.13)
where the errors i are (i.i.d) and normally distributed with zero mean and constant
variance σ2 , and β0 and β1 are unknown regression coefficients which are usually
estimated using least-squares or other similar techniques.
The likelihood function in this case is given by the joint probability of the
75
data:
p(y1, . . . , yn|β0, β1, σ2 ) ∝
[ 1
σ2
]n
2 n∏
i=1
e
[
− 12σ2
(yi−β0−β1ξi)2
]
(4.14)
However, there arises cases where the ξi’s are not directly observed. Instead,
xi’s are observed with a normally distributed additive error νi, which are related by:
xi = ξi + νi ; i = 1, 2, 3, . . . , n. (4.15)
The Equation (4.15) constitutes the classical Error-In-Variable model.
Cheng and Ness [60] states three different sub-classifications of the classical EIV
model based on the assumptions on ξi’s. These are given as follows:
1. Functional Model: ξi’s are assumed unknown constants.
2. Structural Model: ξi’s are (i.i.d) random variables and independent of the error
νi’s. In this case,
E(ξi) = κ; var(ξi) = σ2ξ (4.16)
3. Ultrastructural Model: Similar to the structural model, ξi’s are independent
random variables, but not identically distributed. This means that they have
possibly dissimilar means, κi’s, and a common variance σ2ξ . If κ1 = κ2 = . . . κn,
the ultrastructural model reduces to the structural one, whereas if σ2ξ = 0, the
ultrastructural model reduces to the functional model.
76
For considerations of simplicity the functional model is used for the rest of this
work. To use the data error model together with the EIV model, certain assumptions
are made about the errors i and νi in Eqns. (4.13) and (4.15) respectively. These
are:
E(νi) = E(i) = 0 (4.17a)
var(νi) = σ2ν , var(i) = σ2 ∀i (4.17b)
cov(νi, νj) = cov(i, j) = 0 ∀i 6= j (4.17c)
cov(νi, j) = 0 ∀i, j. (4.17d)
The likelihood function is now the joint distribution of the xi’s from Eqn.
(4.15) and the response variables yi’s from Eqn. (4.13), and is given by:
p(x1, . . . , xn, y1, . . . , yn|β0, β1, σ2 , σ2ν , ξ1, . . . , ξn) ∝
n∏
i=1
[ 1
σ2ν
] 1
2
e
[
− 12σ2ν
(xi−ξi)2
]
×
n∏
i=1
[ 1
σ2
] 1
2
e
[
− 12σ2
(yi−β0−β1ξi)2
]
(4.18)
4.3.3 EIV-based Changepoint Modeling: Likelihood
From Eqn. (4.12), it can be seen that in the non-contact regime, the covariates
ξi = δi, i = 1, 2, . . . , k are fixed and observed; therefore implementation of the EIV
model is not necessary.
77
The joint distribution in the non-contact regime is therefore given as:
p(d1, . . . , dk|β1, σ21, k) ∝
k∏
i=1
[ 1
σ21
] 1
2
e
[
− 12σ21
(di−β11−β12δi)2
]
(4.19)
where the response variables are the deflection data d1, . . . , dk.
In the contact regime, it one can write, using Eqn. (4.12):
ξk+1 = sfk+1
...
ξn = sfn (4.20)
The covariates ξk+1, . . . , ξn in Eqn. (4.20) are multiples of the probe compli-
ance s. When s is observed without errors, the estimation approach proposed in this
work reduces to the regular Bayesian Changepoint model [48],[61], since ξ1, . . . , ξn
are fixed and observed. However, variations in the probe calibration measurements
as evidenced by the results in Table 4.2 indicate that incorporation of some notion
of randomness on the nature of s is necessary.
The variable s is modeled as an unobserved variable which needs to be esti-
mated. Instead of observing s, the random variable sd is observed with an additive
Gaussian error ψ ∼ N(0, σ2sd), which is given by:
sd = s+ ψ; ψ ∼ N(0, σ2sd) (4.21)
78
Using the observed probe compliance sd, one can write:
xk+1 = sdfk+1
...
xn = sdfn (4.22)
where xk+1, . . . , xn are the observed covariates analogous to the observed xi’s in
Eqn. (4.15).
At this point, it would seem natural to define a joint distribution of the data
(xi, di) similar to Eqn. (4.18). However, from Eqn. (4.22) it is clear that xk+1, . . . , xn
are not independent; instead, they are multiples of sd. Likewise, the unobserved
ξk+1, . . . , ξn are all multiples of s. Moreover, it is important to note that the inde-
pendence condition from Eqn. (4.17c) is violated, since:
cov(νi, νj) = cov(xi, xj) = fifjcov(sd, sd)
6= 0 ∀ i 6= j (4.23)
The direct implication of the violation of the independence condition of Eqn.
(4.17c) is that the covariate distribution can simply be expressed in sd and s, instead
of the covariates xk+1, . . . , xn and ξk+1, . . . , ξn. The EIV model in the contact
regime therefore can be reduced to:
p(sd|s, σ2sd) ∝
[ 1
σ2sd
] 1
2
e
[
−(sd−s)
2
2σ2sd
]
(4.24)
79
This leads to considerable reduction in the complexity of the joint distribu-
tion of the contact-regime data. The joint distribution of the EIV model and the
deflection data in the contact regime now can be written as:
p(sd, dk+1, . . . , dn|β2, σ22 , k, s) = p(sd)
n∏
i=k+1
p(di)
∝
[ 1
σ2sd
] 1
2
e
[
−(sd−s)
2
2σ2sd
] n∏
i=k+1
[ 1
σ22
] 1
2
e
[
− 12σ22
(di−β21−sβ22fi)2
]
(4.25)
Combining the joint distribution in the non-contact regime [Eqn. (4.19)] with
the distribution of sd in the contact regime [Eqn (4.25)], the data likelihood for the
whole dataset is given by:
p(sd, d1, d2, . . . , dn|β1,β2, σ21, σ22, k, s) = p(sd)
n∏
i=1
p(di)
∝
[ 1
σ2sd
] 1
2
e
[
−(sd−s)
2
2σ2sd
] k∏
i=1
[ 1
σ21
] 1
2
e
[
− 12σ21
(di−β11−β12δi)2
] n∏
i=k+1
[ 1
σ22
] 1
2
e
[
− 12σ22
(di−β21−sβ22fi)2
]
(4.26)
Based on the candidate contact point k, the deflection data can be re-written in
a compact vector form, as, d1 =
[
d1 d2 . . . dk
]T
, d2 =
[
dk+1 dk+2 . . . dn
]T
and d =
[
d1 d2
]
. The right hand side of [Eqn. (4.12)] form the design matrices
given by:
X1 =


1 1 . . . 1
δ1 δ2 . . . δk


T
80
X2 =


1 1 . . . 1
sfk+1 sfk+2 . . . sfn


T
(4.27)
Using the notations in Eqn. (4.27), the data likelihood of (sd,d) can be re-
written as:
p(sd,d|β1,β2, σ21, σ22, k, s) ∝
[ 1
σ2sd
] 1
2
e
[
− 12σ2sd
(sd−s)2
]
×
[ 1
σ21
] k
2
e
[
− 12σ21
||d1−X1β1||2
]
×
[ 1
σ22
]n−k
2
e
[
− 12σ22
||d2−X2β2||2
]
(4.28)
4.3.4 EIV-based Changepoint Modeling: Posterior
Table 4.3: Conjugate Prior Distributions
Model Distribution Hyperparameters
Parameter Family
k Uniform ∼ U(1, n)
β1,β2 Normal ∼ N(β¯1, Λ¯
−1
1 ) , ∼ N(β¯2, Λ¯
−1
2 )
σ21 , σ22 Inverse Gamma ∼ IG(a0, b0) , ∼ IG(a0, b0)
s Normal ∼ N(µsp, σ2sp)
The posterior distribution, expressed as posterior distribution ∝ data likelihood×
prior distribution, and can be written as:
p(β1,β2, σ21, σ22 , k, s|sd,d, σ2sd) ∝p(d, sd|β1,β2, σ
2
1, σ22, k, s)× pi(β1)× pi(β2)×
pi(σ21)× pi(σ22)× pi(k)× pi(s) (4.29)
where pi(.) indicates the prior distribution of the variable in parenthesis. Conjugate
81
priors (Table 4.3) are used to ensure that the posterior distribution is separable [61].
Combining the distribution of the data with the priors, the posterior can be written
as:
p(β1,β2, σ21, σ22, k, s|sd,d, σ2sd) ∝
[ 1
σ21
] k
2
[ 1
σ22
]n−k
2
[ 1
σ2sd
] 1
2
e
[
− 12σ21
||d1−X1β1||2
]
e
[
− 12σ22
||d2−X2β2||2
]
e
[
− 12σ2sd
(sd−s)2
]
e
[
− 12 (β1−β¯1)T Λ¯1(β1−β¯1)
]
e
[
− 12 (β2−β¯2)T Λ¯2(β2−β¯2)
]
(σ21)−a0−1e
[
− b0
σ21
]
(σ22)−a0−1e
[
− b0
σ22
][ 1
σ2sp
] 1
2
e
[
− 12σ2sp
(s−µsp )2
]
(4.30)
A schematic of the EIV-Bayesian Changepoint algorithm is given in Fig. 4.4.
Data Likelihood
3 3.5 4 4.5
−20
0
20
40
Raw AFM data
Probe Calibration 
data
Priors
Joint Posterior GibbsSampling Marginal
Posteriors
Force-Indentation 
relationship
Figure 4.4: Schematic of the EIV-Bayesian Changepoint algorithm c©2014 IEEE.
4.3.5 EIV-based Changepoint Modeling: Gibbs Sampling
The marginal posteriors are obtained from the joint posterior distribution of
Eqn. (4.30) using Gibbs Sampling, which is a Markov Chain Monte Carlo (MCMC)
sampling technique wherein the marginal posterior in each parameter of the joint
82
posterior is derived by separating out the terms in the joint posterior corresponding
to each parameter and conditional upon the rest. The posterior is sampled iteratively
for the parameters (s,β1,β2, σ21, σ2s , k).
While sampling for s(i+1) given i previous iteratively sampled steps, it can be
observed that s appears in the data likelihood and the prior for s. The terms in the
posterior distribution containing s are given by:
p(s) ∝ e
[
− 12σ22
||d2−X2β2||2
]
e
[
− 1
2σ2sd
(sd−s)2
]
e
[
− 12σ2sp
(s−µsp )2
]
∝ e
[
− 12σ22
n∑
j=k+1
{di−β21−β22sfj}2
]
e
[
− 12σ2sd
(sd−s)2
]
e
[
− 1
2σ2sp
(s−µsp)2
]
(4.31a)
The coefficients of s2 and s can be combined to complete the square in s by a simple
rearrangement of terms in Eqn. (4.31a). Since the constant term in the square in s
is not dependent on s, it can be ignored without altering the Gibbs Sampling step.
Therefore:
s(i+1) ∼ p(s|β(i)1 ,β
(i)
2 , σ21
(i), σ22
(i), k(i))
∼ N(µ˜s, σ˜2s)
83
where :
µ˜s =
[(σ2spsd + σ2sdµsp
σ2sp + σ2sd
) σ22
β222
n∑
j=k+1
f 2j
+
(
n∑
j=k+1
(dj − β21)fj
β22
n∑
j=k+1
f 2j
)( σ2spσ2sd
σ2sp + σ2sd
)][
σ22
β222
n∑
j=k+1
f 2j
+
( σ2spσ2sd
σ2sp + σ2sd
)]−1
(4.31b)
σ˜2s =
[
β222
σ22
n∑
j=k+1
f 2j +
1
σ2sd
+ 1σ2sp
]−1
(4.31c)
Likewise, the regression parameters β1 and β2 appear in multivariate normal form in
the data likelihood and the priors for β1 and β2 respectively, and the same strategy
is employed to complete the square and factor it out of the posterior. The terms in
the posterior distribution containing the regression parameters are given by:
p(βm) ∝ e
[
− 12σ2m
||dm−Xmβm||2
]
e
[
− 12 (βm−β¯m)T Λ¯m(βm−β¯m)
]
∝ e
[
− 12σ2m
{||dm−Xmβˆm||2+(βm−βˆm)TXTmXm(βm−βˆm)}
]
e
[
− 12 (βm−β¯m)
T Λ¯m(βm−β¯m)
]
(4.32a)
where m = 1, 2 for the non-contact and the contact regimes respectively and βˆm =
(XTmXm)−1XTmdm is the least-squares solution.
The following identity is used to complete the square in the multivariate case
84
[62]:
2∑
p=1
(α−αp)TΣp(α−αp) = (α− α˜)TΓ(α− α˜) + f(αp,Σp) (4.32b)
whereΣ1,Σ2 are positive semi-definite symmetric matrices, α˜ = (Σ1+Σ2)−1(Σ1α1+
Σ2α2) and Γ = Σ1 +Σ2.
Using Eqn. (4.32b), one arrives at the following expression:
β(i+1)1 ∼ p(β1|s(i+1),β
(i)
2 , σ21
(i), σ22
(i), k(i))
∼ N(β˜1, Σ˜1) (4.32c)
β(i+1)2 ∼ p(β2|s(i+1),β
(i+1)
1 , σ21
(i), σ22
(i), k(i))
∼ N(β˜2, Σ˜2) (4.32d)
where
β˜m = (σ−2m XmTXm + Λ¯m)
−1(σ−2m XmTdm + Λ¯mβ¯m) (4.32e)
Σ˜m = (σ−2m XmTXm + Λ¯m)
−1 (4.32f)
for m = 1, 2.
The variances σ21 and σ22 appear in the data likelihood and the priors for the
variances. Likewise, the posterior can be factorized to contain terms only in σ21 and
85
σ22 and the following expression results:
p(σ21) ∝
[ 1
σ21
]k
2
e
[
− 12σ21
||d1−X1β1||2
][ 1
σ21
]a0+1
e
[
− b0
σ21
]
∝
[ 1
σ21
]k
2+a0+1
e
[
− 1
σ21
(
0.5||d1−X1βm||2+b0
)]
(4.33a)
and
p(σ22) ∝
[ 1
σ22
]n−k
2
e
[
− 12σ22
||d2−X2β2||2
][ 1
σ22
]a0+1
e
[
− b0
σ22
]
∝
[ 1
σ22
]n−k
2 +a0+1
e
[
− 1
σ22
(
0.5||d2−X2β2||2+b0
)]
(4.33b)
Hence, one can write:
σ2,(i+1)1 ∼ p(σ21|s(i+1),β
(i+1)
1 ,β
(i+1)
2 , σ22
(i), k(i))
∼ IG(a˜1, b˜1) (4.33c)
σ2,(i+1)2 ∼ p(σ22|s(i+1),β
(i+1)
1 ,β
(i+1)
2 , σ
2,(i+1)
1 , k(i))
∼ IG(a˜2, b˜2) (4.33d)
where
a˜1 = a0 + k/2 (4.33e)
a˜2 = a0 + (n− k)/2 (4.33f)
b˜m = b0 + 0.5(||dm −Xmβm||)2 (4.33g)
86
for m = 1, 2.
In the case of the contact point index, a flat uninformative prior [pi(k) = 1]
is used to avoid biasing the posterior. The posterior in k is therefore the discrete
distribution function obtained by varying k over the data likelihood, and is given
by:
k(i+1) ∝ p(k|s(i+1),β(i+1)1 ,β
(i+1)
2 , σ
2,(i+1)
1 , σ
2,(i+1)
2 )
∝
[ 1
σ21
]k
2
e
[
− 12σ21
||d1−X1β1||2
][ 1
σ22
]n−k
2
e
[
− 12σ22
||d2−X2β2||2
]
(4.34)
4.4 Results and Discussions
4.4.1 Implementation of the Gibbs Sampler
Excepting the contact point index k, all the other marginal posteriors can be
directly sampled from their respective family of parametric probability distributions
due to the use of conjugate priors. Rejection sampling [63] is used to sample for k
[Eqn. (4.34)] to complete the Gibbs Sampling step.
The initial 3000 iteration results of the Gibbs Sampler are rejected, which
constitutes the burn-in period. Sampling is terminated when increasing the number
of iterations do not alter the nature of the marginal posteriors. Typically, sampling
is terminated after 200,000 iterations. With respect to the probe compliance, sd is
set to a pre-experiment calibration value (= 16.1 m/N for Probe 1 and 0.339 m/N
for Probe 2). The hyperparameter values of the probe compliance are set to those
obtained in Table 4.2 , i.e. µsp = 16.38 m/N, σ2sp = 1.632 m2/N2 for Probe 1 and
87
Table 4.4: Hyperparameter Values
Hyperparameters Values
β¯1, β¯2
[
0 0
]T
Λ¯1, Λ¯2
[
10−15 0
0 10−15
]
(a0, b0) (1, 2) (Hertz Model)
(1, 0.002) (Long’s Model)
(µsp, σ2sp) (16.38, 1.632) (for Probe 1)
(0.328, 0.0322) (for Probe 2)
µsp = 0.328 m/N, σ2sp = 0.0322 m2/N2 for Probe 2. Using such an informative prior
based on Table 4.2 makes it possible to incorporate prior knowledge of multiple
probe calibrations. The prior and data variances of the probe compliance are also
assumed to be equal i.e. σ2sp = σ2sd . This allows us to give equal weightage to
the prior and data values. Setting the data and prior variances equal also avoids
estimation problems, which are typical when the data variance is unknown [64].
With respect to the other parameters, no useful information is available beforehand,
consequently, uninformative priors are assigned to them. The contact point k is
assigned a uniform distribution, U(1, n), while β1,β2, σ21 and σ22 are assigned disperse
prior distributions. The hyperparameter values used during sampling is given in
Table 4.4.
4.4.2 Simulated Data
Owing to lack of robust experimental methods to determine the contact point
from AFM force curves, the EIV-Bayesian Changepoint algorithm is first validated
88
0 100 200 300 400 500 600 700 800 900 1000
30
35
40
45
50
55
60
65
70
630 635 640 645 650 655 660 665 670 675 680
0
0 .02
0 .04
0 .06
0 .08
0 .1
0 .12
0 .14
0 .16
0 .18
Posterior in  k
%
 F
re
qu
en
cy
Y
Contact Point Index 
Contact Point Index 
Posterior 
Mean
Ground Truth
95% Confidence
Interval 
Ground
 Truth
(a) 
(b) 
Figure 4.5: Results of the EIV-Bayesian Changepoint algorithm on simulated data
with n = 1000 datapoints. (a) Simulated AFM force curve with contact occurring
at datapoint k = 650 and (b) marginal posterior of the contact point index k, along
with the posterior mean bkˆc = 653 and the 95% confidence intervals [642, 660].
89
on a simulated AFM force curve with a known contact point. The Hertz model (with
power law 1.5) is used as the force-indentation relationship, and the calibration re-
sults for Probe 1 are used to specify the probe compliance data and hyperparameter
values. The following parameters were used to generate the entire data:
• Datapoints n = 1000. Contact Point Index k = 650.
• Pre-contact regime standard deviation σ1 = 1.20. Contact regime standard
deviation σ2 = 1.50.
• Pre-contact regime regression parameters β1 =
[
β11 β12
]T
=
[
35 0.001
]T
.
• Contact regime regression parameters β2 =
[
βj1 βj2
]T
=
[
35.01 0.3
]T
.
• Observed probe compliance sd = 16.10. Observed probe compliance variance
σ2sd = 1.632.
The simulated AFM curve is shown in Fig. 4.5(a). The noise parameters
σ1 and σ2, were selected to mimic experimental AFM force curves obtained on
extremely compliant live cells, where the contact uncertainty is more severe than
other biological specimens.
The result of sampling the marginal posterior in the contact point index k
according to Eqn. (4.34) is shown in Fig. 4.5(b). The true contact point is at
datapoint 650, while the posterior mean of the contact point is given by bkˆc = 653,
a very negligible error given the noisy nature of the transition. Moreover, it is sig-
nificant to note that the 95% confidence interval (C.I.) given by [642, 660], contains
the true contact point. This is one of the advantages of using an interval based
90
0.25 0.3 0.35
0
0.02
0.04
0.06
0.08
0.1
0.12
14 16 18 20
0
0.02
0.04
0.06
0.08
0.1
0.12
1.3 1.4 1.5 1.6 1.7 1.8
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
1.1 1.2 1.3 1.4 1.5
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
%
 F
re
qu
en
cy
%
 F
re
qu
en
cy
%
 F
re
qu
en
cy
%
 F
re
qu
en
cy
Probe Compliance 
Pre-contact Standand Deviation Contact Standard Deviation
Ground
Truth
(a) (b)
(c) (d)
Posterior
Mean
Posterior
Mean
Ground
Truth
Posterior
Mean
Ground
Truth
Posterior
Mean
Figure 4.6: Marginal posterior distributions in (a) β22 (b) s (c) σ1 and (d) σ2.
91
Table 4.5: Model performance on simulated data
Ground Truth Posterior Means
k 650 653
β22 0.3 0.306
σ1 1.20 1.24
σ2 1.50 1.50
(sd, µsp) (16.10,16.38) 16.74 (sˆ)
method over point based methods in estimating the contact-point and the mechan-
ical properties: even though the true contact point was missed by 3 datapoints, the
95% C.I contained the true contact point. It is also evident from Figs. 4.5(a) and
4.5(b) that integration of the EIV model into a standard Changepoint model does
not adversely affect Changepoint model’s ability to estimate the contact point.
The marginal posteriors in s, β22, σ21, σ22 are given in Fig. 4.6, and the posterior
means are given in Table 4.5. From the results, it can be seen that the EIV-Bayesian
Changepoint model performs quite satisfactorily on simulated data. The estimation
error in the regression coefficient of interest, β22 is 2%, which is quite reasonable
given the noisy nature of the simulated AFM force curve. The posterior means of
σ1 and σ2 are also within 3.5 % of the ground truth values. The probe compliance
s shows a posterior mean (sˆ = 16.74) with a dispersed distribution due to the
variability introduced by the variance terms σ2sd and σ2sp.
4.4.3 AFM studies on histological breast tissue
This section deals with the implementation of the algorithm on AFM indenta-
tion datasets obtained from histological breast tissue samples. The AFM indenta-
92
Tracking
ROI
Probing
 ROI
AFM
Probe
Figure 4.7: Image of tissue region where probing was conducted to obtain the results
in this section. AFM indentation experiments on breast tissue were conducted in
the probing ROI c©2014 IEEE.
tion data was generated from raster scanning of the probing ROI shown in Fig. 4.7.
Probe 2 was used for the AFM experiments, and Long’s contact model [59] was used
as force-indentation relationship.
Due to the absence of existing algorithms that quantify the elastic modulus
estimation errors due to both contact uncertainty and probe calibration variability,
it is difficult to compare the efficacy of the algorithm in reflecting the variability in
the estimated elastic modulus due to these two sources of uncertainty. To overcome
this, the following steps are implemented:
1. The 95% bounds from multiple calibration experiments in Table 4.2, i.e.
s¯2.5% = 0.328− 2(0.032) m/N and s¯97.5% = 0.328 + 2(0.032) m/N are used to
set deterministic bounds of the probe compliance variability from the probe
compliance.
2. The posterior mean of the contact point, kˆ, is obtained through a standard
93
Bayesian Changepoint approach [48] (without any variations in the probe com-
pliance).
3. Using s¯2.5% and s¯97.5% and the contact point posterior mean from step (2), the
bounds of the elastic modulus, E¯2.5% and E¯97.5%, are estimated using least-
squares fitting to the Long’s model. These bounds are then compared to
the 95% confidence intervals of the marginal posterior in the elastic modulus
obtained from the EIV-Bayesian Changepoint approach.
The results of the EIV-Bayesian Changepoint algorithm applied on a repre-
sentative AFM force curve on tissue specimens are presented in Figs. 4.8 and 4.9.
The AFM deflection data together with the posterior mean of the contact point
(bkˆc = 433) is shown in Fig. 4.8(a), and the marginal posterior in k is shown in
Fig. 4.8(b).
The marginal posterior in s is shown in Fig. 4.9(a). The marginal posterior
in β22 [Fig. 4.9(b)] is easily transformed into the elastic modulus, E, due to the
linear relationship between them [Eqn. (4.10)]. The posterior mean of the elastic
modulus, Eˆ = 119.01 kPa compares favorably with previously reported values on
histological breast tissue [65].
Also displayed in the Figs. 4.9(a) and (b) are the bounds (s¯2.5%, s¯97.5%) and
(E¯2.5%, E¯97.5%) respectively. It is evident that the 95% confidence intervals of the
elastic modulus are well bounded by E¯2.5% and E¯97.5%. Indeed, this is one of the
advantages of using an integrated EIV-Bayesian Changepoint approach: lower vari-
ability in the estimated elastic modulus compared to using bounds on the calibration
94
0 200 400 600 800
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7(b)
%
 F
re
qu
en
cy
Contact point Index 
0 200 400 600 800
−10
0
10
20
30
40
50
60
Posterior
Mean 
(a)
De
fle
ct
io
n 
(n
m
)
Datapoints
420 430 440 450
−0.55
−0.5
−0.45
−0.4
−0.35
−0.3
Figure 4.8: Results of the EIV-Bayesian Changepoint algorithm on a representative
AFM force curve on breast tissue with n = 783 datapoints. (a) AFM force curve
with the posterior mean of the contact point index at bkˆc = 433. (b) Marginal
posterior of k.
95
0.25 0.3 0.35 0.4 0.45
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
90 100 110 120 130 140 150 160
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
Elastic Modulus      (kPa)Probe Compliance      (m/N)
%
 F
re
qu
en
cy
%
 F
re
qu
en
cy
(a) (b)
Posterior Mean 
95% Confidence
 Interval
Posterior Mean 
95% Confidence 
Interval
Figure 4.9: (a) Marginal posterior in s and (b) marginal Posterior of E.
data and point estimates of the contact point to compute the elastic modulus.
4.4.4 Sensitivity Analysis
The dispersed nature of the marginal posterior in the elastic modulus in
Fig. 4.9(b) raises an interesting question - what are the individual contributions of
the contact point uncertainty and the probe compliance variability in the marginal
posterior distribution in E? More importantly, does the marginal posterior in E
reflect changes in the hyperparameter values of the probe compliance, given that an
informative prior was used for it?
4.4.4.1 Varying σsp
The effect of using hypothetical values of σsp on the marginal posteriors is
illustrated in Fig. 4.10. The dispersion of the marginal posterior in E increases
96
0.000 0.012 0.022 0.032 0.042 0.052
80
100
120
140
160
0.000 0.012 0.022 0.032 0.042 0.052
0.25
0.3
0.35
0.4
0.45
0.000 0.012 0.022 0.032 0.042 0.052
426
428
430
432
434
436
 (k
Pa
)
(m
/N
)
(a)
(b)
(c)
(m/N)
97.5%
2.5%
mean
max
min
Figure 4.10: Boxplots of the marginal posterior distributions (a) in E, (b) in s and
(c) k. Please note that the σsd has been assumed to be equal to σsp. The fourth
boxplot (σsp = 0.032) corresponds to the experimental probe compliance variance
as obtained in Table 4.2 c©2014 IEEE.
97
with σsp , as shown in Fig. 4.10(a). This is an expected observation, since larger
variations in the calibrated probe compliance would lead to greater uncertainty in the
elastic modulus of the probed regions. The individual contribution of contact point
uncertainty on marginal posterior in E is given in the first boxplot of Fig. 4.10(a)
(corresponding to σsp = 0.000).
Likewise, increase in σsp leads to greater dispersion in the marginal posteri-
ors in s [Fig. 4.10(b)]. It is worth noting that changes in σsp does not cause any
substantial change to the posterior means of E or s.
The boxplots of Fig. 4.10(c) shows the effect of σsp on the marginal posterior
in k. It can be observed that changes in σsp has no significant effect on the posterior
in k - the 95% confidence intervals varies between the 428th and 435th datapoints.
This is an important observation, since uncertainties in the probe compliance and
the contact point are physically unrelated - probe compliance variations are the
result of calibration errors, while contact point uncertainty occurs because of the
soft nature of biological tissues and the absence of a perceptible attractive region in
the AFM force curve in liquid.
4.4.4.2 Varying µsp
The effect of varying the hyperparameter µsp on the marginal posteriors in
E and s are shown in Figs. 4.11(a) and 4.11(b) respectively. The previous values
of sd(= 0.339) m/N and σsd(= 0.032) m/N and the assumption that σ2sd = σ2sp
are retained. The first boxplot shows the marginal posteriors when µsp is set to
98
0.222 m/N, the probe manufacturer’s rated compliance. This is often useful when
the AFM user performs a single probe calibration experiment, and wishes to see the
result of using the rated compliance instead of repeating the calibration experiments.
The second boxplot uses µsp = 0.328 m/N, obtained from Table 4.2.
Reducing the hyperparameter µsp from 0.328 m/N to 0.222 m/N leads to an
increasing trend in the marginal posterior in E [Fig. 4.11(a)] and a corresponding
decreasing trend in the marginal posterior in s [Fig. 4.11(b)]. This is understand-
able, given the inverse relationship between s and E [see Eqn. (4.10)]. This inverse
relationship is also responsible for the increased dispersion in the first boxplot of
Fig. 4.11(a) compared to the second; in contrast both boxplots in Fig. 4.11(b) show
largely similar dispersed behavior.
It is also evident from Fig. 4.11(b) that the probe compliance posterior means
(sˆ) lie approximately midway between the compliance data sd and the hyperparam-
eter µsp (sˆ = 0.283 m/N for µsp = 0.222 m/N and sˆ = 0.333 m/N for µsp = 0.328
m/N). This is a direct consequence of using equal weightage to the data and prior
variances of the probe compliance, i.e. (σ2sd = σ2sp).
4.5 Conclusions
This chapter outlines the formulation of an integrated EIV based-Bayesian
Changepoint algorithm used to obtain robust estimates of the elastic modulus of
biological specimens from AFM force curves in the presence of variability in the
contact point and probe spring constant. The major highlights of this chapter are
99
0.222
0.2
0.25
0.3
0.35
0.4
(m/N)
0.222
100
120
140
160
0.328
(a)
(b)
0.328
(k
Pa
)
(m
/N
)
97.5%
2.5%
mean
max
min
Figure 4.11: Boxplots indicating the effect of varying µsp on the marginal posteriors
(a) in E and (b) in s c©2014 IEEE.
enumerated as follows:
1. Transforming the probe spring constant to the equivalent probe compliance,
a standard Bayesian Changepoint problem [48] could be set up to accommo-
date variability in the probe spring constant. This allowed the application
of error-in-variables based methods to be applied to the Bayesian Change-
point algorithm. Multiple probe calibration experiments were conducted to
obtain the variability in the probe compliance using the thermal calibration
method. These results were used to specify the hyperparameters of the prior
distribution of the probe compliance.
2. The proposed algorithm performed satisfactorily on simulated AFM data, and
the estimation errors were within 4%. Also, the ground truth was contained
within 95% confidence intervals of the marginal posteriors.
100
3. On implementing the algorithm on AFM datasets obtained from breast tissue
specimens, it was observed that using an integrated approach, the variability
in the estimates of the elastic modulus could be reduced, as opposed to point-
based approaches to estimate the elastic modulus by treating each source of
variability piecemeal.
4. A sensitivity analysis on the parameters of the posterior distribution showed
that the proposed algorithm model responded satisfactorily to a wide range of
hyperparameter values in the probe compliance.
Accurate estimation of the contact point shown in Fig. 4.5 and the sensitivity
plots shown in Figs. 4.10 and 4.11 confirm the research hypothesis of this chapter
that EIV-Changepoint algorithm successfully estimates the elastic modulus due to
the aforementioned sources of uncertainty.
One of the appealing aspects of using the EIV-Bayesian Changepoint formula-
tion is that it makes it possible to automate data analysis of AFM force curves, since
the contact point from the AFM force curves can be determined fairly accurately
without subjective user input. As mentioned previously in chapters 2 and 3, AFM
indentation studies on breast tissue for histopathological inference necessitates very
large scale AFM indentation data acquisition and processing. A robust and com-
putationally efficient approach can significantly improve throughput of AFM based
characterization of biomaterials by allowing batch-processing of AFM curves.
The proposed approach also obviates the need to truncate AFM force curve
datasets to ensure applicability of linear elastic theory. The use of Long’s contact
101
model [59] accounts for geometric and material nonlinearities and therefore the whole
AFM force curve can be processed to estimate material properties.
However there are certain limitations of this work that require elaboration.
These are stated as follows:
1. The EIV-Bayesian Changepoint algorithm assumes that the force-indentation
contact relationship appropriately fits the post-contact AFM deflection data.
As a result, the algorithm has been implemented only on those AFM datasets
where the post-contact deflection data could be well described by the contact
model in question, i.e., Hertz or Long’s contact model. For force curves that
display sharp nonlinearity, Long’s contact model is not likely to be an ap-
propriate contact model since it only accounts for neo-Hookean hyperelastic
behavior - which is a first order approximation of higher order hyperelastic
material models. The effect of contact estimation on force curves with pro-
nounced nonlinear effects are dealt with separately in Chapter 5.
2. Errors in the observed probe compliance have been assumed to follow a Gaus-
sian distribution - this assumption has not been experimentally validated. The
assumption of normality in the observed probe compliance and in the prior
distributions of the probe compliance was primarily aimed at simplifying the
sampling of the joint posterior using standard MCMC methods like Gibbs
Sampling. Indeed, large scale probe calibration experiments are necessary to
ascertain the nature of the distribution to which the calibrated probe compli-
ance errors belong.
102
3. In the two regime regression model (Eqn. 4.12) used for the EIV-Bayesian
Changepoint analysis, continuity was not assumed at the contact point. This
did not impact the performance of the EIV-Bayesian Changepoint algorithm,
primarily because the algorithm was only implemented with force curves that
could be described well by the contact model chosen in the analysis (Hertz
or Long’s [59] contact model). For force curves that cannot be adequately
described by the contact model, the contact point estimates are likely to be
erroneous. These issues are dealt with in detail in chapter 5.
103
Chapter 5
Constitutive material modeling of the tissue specimens under AFM
indentation
5.1 Overview
In chapters 2 and 3, the Hertzian contact model was used to extract the tis-
sue elastic modulus, which relies on the theory of linear elasticity and infinitesimal
strains, to estimate the elastic modulus of the tissue specimens. This was on account
of its computational ease and also because the primary objective was to investigate
any elasticity differences between normal and cancerous specimens, not so much to
assess the appropriateness of the contact model for any given force curve. In chap-
ter 4, Long’s contact model [59] was implemented in the EIV-Bayesian Changepoint
estimation algorithm, which relaxes the assumptions of linear elasticity in the ma-
terial by accounting for neo-Hookean hyperelasticity. However, the algorithm was
only implemented on those force curves that could be adequately described by Long’s
contact relationship.
In this chapter, greater emphasis is placed on the choice of the contact model to
describe AFM force curves with varying levels on nonlinear elasticity and numerical
approaches are explored that might be better suited for elastic modulus estimation
from the tissue specimens. The coupled problem of contact estimation from highly
104
nonlinear force curves is also investigated, and a modification to the contact estima-
tion algorithm from chapter 4 is proposed, which can be used in conjunction with
numerical techniques for tissue mechanical property estimation.
It should be mentioned here that in this chapter, the probe spring constant (or
compliance) is assumed to be observed without errors and no consider variability in
the probe’s spring constant is considered in the work outlined in this chapter.
5.2 Experimental Challenges in enforcing Hertzian assumptions
The applicability of the Hertz contact model for analyzing AFM indentation
data on biomaterials has critically studied in the past [1],[66] - nevertheless, a critical
assessment of each of the assumptions in Hertzian theory in the context of prevail-
ing experimental conditions during AFM indentation on breast tissue specimens is
presented in the following discussion.
It may be recalled from chapter 2 that the assumptions in Hertzian contact
model for spherical indentation of tissue specimen are:
(H1) The contact between the indenter and the specimen is considered frictionless.
(H2) There is no adhesion between the sphere and the indenter.
(H3) The displacements and strains in the specimen are assumed to be small so
that a linearized continuum theory can be used.
(H4) The specimen is assumed to be an infinite half-space.
(H5) The specimen is linearly elastic.
105
(H6) The specimen is homogeneous and isotropic.
The assumption (H1) of frictionless contact implies that shear tractions at the
contact interface are zero. Since AFM indentation on tissue specimens are primarily
conducted in PBS medium, this is a reasonable assumption to make in the present
study.
11.5 12 12.5 13 13.5 14 14.5
−300
−250
−200
−150
−100
−50
0
Z-Position,  
De
fle
ct
io
n,
 
12.8 12.9 13 13.1 13.2 13.3 13.4
−300
−295
−290
−285
−280
−275
−270
−265
−260
DisplacementDisplacement
Fo
rc
e
Fo
rc
e
(a) (b)
(c)
No-adhesion Adhesion
Representative AFM force curve
Figure 5.1: Schematic of force-displacement behavior for (a) adhesion-free inden-
tation, (b) indentation with adhesion and (c) representative AFM force curve on
breast tissue specimens showing absence of prominent adhesive region.
The assumption (H2) of adhesion between the tissue and the indenter is not a
trivial one. Live mouse fibroblasts [67] and polyvinyl alcohol gels [68] have previous
shown adhesive behavior when subjected to AFM indentation, which causes the
AFM probe to deflect toward the sample surface resulting in attractive (negative
force) regime in the AFM force curve [see Fig. 5.1(a) and Fig. 5.1(b)]. However, as
106
seen in chapters 2, 3 and 4, no adhesion was observed in the force curves on the
tissue specimens [Fig. 5.1(c)], consequently, modeling adhesion is not a significant
concern in this study.
Figure 5.2: Finite-indentation of a thin specimen by a rigid hemispherical indenter.
Assumptions (H3) and (H4) constitute two of the most critical assumptions
that are often violated during AFM indentation studies on tissue specimens. The
infinitesimal assumptions in (H3) and (H4) imply the following:
∆ R (5.1a)
∆ h (5.1b)
R h (5.1c)
h→∞ (5.1d)
In order to ensure applicability of Eqns. (5.1a) and (5.1b), AFM force curves need to
be acquired in indentation-controlled mode to restrict the sample indentation to the
Hertzian regime. However, conducting AFM experiments in indentation-controlled
107
mode is significantly harder compared to force-controlled mode due to the challenges
in online contact estimation, as stated previously in chapter 2. Coupled with the
need to perform raster-scans on the tissue specimens for large-scale AFM force
curve acquisition, indentation-controlled mode becomes a far less attractive mode
compared to force-controlled mode for high throughput purposes. The drawback
associated with force-controlled AFM experiments, however, is that highly compliant
regions in the tissue register very large indentations, which violate assumptions
(5.1a) and (5.1b).
Applicability of Eqns. (5.1c and 5.1d) can be achieved by adjusting the bead
size with respect to the specimen thickness. However, as observed previously in [24],
probe tips with low radius (for example, pyramidal tips with nm-range radius) can
cause permanent damage to the tissue samples [69]. Even worse, sharp tips possess
the undesirable tendency of including contributions from the substrate [70], leading
to inaccurate material characterization. An optimum range of R ∼ 5−10µm ensures
low stress concentration in the tissue specimens, and also aids on averaging out local
tissue heterogeneity. Tissue specimen thickness also can be varied from h = 4µm
to h = 10µm, however, thick tissue specimens also lead to loss of continuity in the
tissue architecture across consecutive tissue slices. Even in the most optimistic case
of (R, h) = (5, 10) µm, R/h = 0.5, which is clearly too large to enforce infinitesimal
assumptions.
As observed in chapters 2 and 3, the tissue specimens studied in this dis-
sertation consists of complex three-dimensional architecture with functional breast
epithelial cells that are enclosed by collagen-rich stromal tissue - which violates
108
isotropic and homogeneity assumptions [Assumption (H6)]. Moreover, macroscale
breast tissue specimens have previously shown nonlinear elasticity [Assumption
(H5)] [71].
The combined effect of violations in assumptions (H3-H6) is well-illustrated in
Fig. 5.3, which illustrates the spatial distribution of the tissue stiffness [Fig. 5.3(a)]
and the corresponding goodness of fit [Fig. 5.3(b)] on a tissue ROI. It can be observed
that while the tissue stiffness shows clear delineation from the upper left to the
lower right end of the elasticity map, no such patterns can be observed in the
R2 distribution. Figs. 5.3(c) and 5.3(d) show the force-indentation plots of two
neighboring sampled points with similar stiffness, but with vastly dissimilar nature
of the Hertzian fit. These results underscore the need to have specific material
constitutive laws to extract accurate elasticity measures from AFM indentation
data. In the following section, prior work in the development of contact models that
overcome the limitations of the Hertz contact model are critically reviewed.
5.3 Related Work
Previous research addressing non-Hertzian behavior in biomaterials undergo-
ing AFM indentation can be classified by the individual Hertzian assumptions (H1-
H6) that are targeted and corrected. Since the scope of this research is primarily
directed towards addressing violations in assumptions (H3-H6), a description of
some of the approaches proposed in literature that focus on these assumptions only
is presented in the following discussion. Readers are requested to look elsewhere for
109
0.9
0.92
0.94
0.96
0.98
1
100
200
300
400
500
600
700
800
900
1000
1100
1200
Indentation (nm)
In
de
nt
at
io
n 
(n
m
)
Fo
rc
e 
(n
N)
Fo
rc
e 
(n
N)
Indentation (nm)
0 100 200 300 400
0
50
100
150
200
0 100 200 300
0
50
100
150
200
Experiment
Hertz Fit
Experiment
Hertz Fit
(a)
(b)
(c) (d)
kPa
Figure 5.3: Effect of heterogeneity in the spatial stiffness distribution. (a) An elas-
ticity map in the cancer core [from Fig. 3.11(b)], (b) corresponding 2-D representa-
tion of the R-square of Hertzian fits and (c,d) representative force-indentation plots
showing the inadequacy of the Hertz model to describe the entire ROI’s constituent
behavior.
110
corrections to violations in (H1) [72] and (H2) [67].
• Finite Indentation (H3,H4): Improvements to the Hertz contact model to
account for thin specimen indentation was primarily motivated by AFM in-
dentation studies on thin polymer gels, wherein it was observed that the Hertz
contact model significantly overestimated the elastic modulus due to the effect
of the underlying stiff substrate [70]. The most seminal contribution in this
regard points to the work of Dimitriadis and his co-workers [44]. Dimitriadis
et al. [44] appended correction terms to the Hertz contact model that led to
an analytical force-indentation relationship as given below1:
F = 4E
√
R∆1.5
3(1− ν2)
[
1−2α0pi κ+
4α20
pi2 κ
2− 8pi3
(
α30+
4pi2
15 β0
)
κ3+16α0pi4
(
α30+
3pi2
5 β0
)
κ4
]
(5.2)
where E,R,∆, ν have their usual meanings, κ = R∆/h2 and α0, β0 are con-
stants that depend on the Poisson’s ratio ν.
In the limit κ → 0, Eqn. (5.2) reduces to the classical Hertz contact model.
It is noteworthy here that for this to happen, it is not sufficient to restrict the
indentation depth ∆ for a given specimen thickness h only (presumably using
indentation-controlled AFM indentation experiments). Even for moderately
low ∆ values, a large bead radius R could possible make κ non-zero, thereby
invalidating Hertzian assumptions of infinitesimal strains. In the recent past,
Santos et al. [73] have investigated the influence of probe geometry (with
pyramidal/conical indenters of nm-range radius) on the elastic modulus from
1Validation of the Dimitriadis model using finite element analysis is given in Appendix A.
111
thin specimens using numerical methods.
• Nonlinear Elasticity (H5): To account for nonlinear elasticity effects in bioma-
terials, Costa et al. [1] proposed a “pointwise modulus” approach to estimate
the elastic modulus as a function of the indentation in the sample. For a
purely Hertzian response, the “pointwise modulus” is constant during the en-
tire indentation regime, whereas for violations in Hertzian assumptions, the
“pointwise modulus” show distinct nonlinearity. While this approach success-
fully differentiates a Hertzian response from a non-Hertzian one, the effect
of material and geometric nonlinearities (finite indentation) are coupled and
there is no way to isolate these two violations in Hertzian assumptions using
the “pointwise modulus” approach solely.
Drawing upon parallels between uniaxial compression and indentation, Lin
et al. [74][75] proposed a reduction formulation to convert force-indentation
datasets into “indentation stress” (σ∗) v/s “indentation strain” (∗) plots. This
was accomplished by dividing the contact force F by the contact area to obtain
σ∗, and expressing ∗ as a function of the contact area [76], as shown below
[74]:
σ∗ = Fpia2 (5.3a)
∗ = 0.2 aR (5.3b)
a = R1/2∆1/2 (5.3c)
112
where a is the contact area. For a purely Hertzian response, σ∗ - ∗ are given
by [74]:
σ∗ = 20E3pi(1− ν2)
∗ (5.4)
It is noteworthy that σ∗ varies linearly with ∗ for Hertzian indentation (Eqn.
5.4), which is expected due to assumptions of material linearity in Hertzian
theory. Lin [74] then demonstrated that force-indentation equations could
be obtained for a wide range of hyperelastic materials by obtaining σ∗ as a
function of ∗ from their strain energy densities and then recasting them into
Eqns. (5.3a) and (5.3b).
While this approach is attractive due to its mathematical elegance in defining
explicit force-indentation relations for hyperelastic materials, its fundamental
shortcoming is that finite-indentation effects are not accounted for. This is
well demonstrated in Fig. 5.4. In Fig. 5.4(a), the purple force-indentation
plot corresponds to a simulated Hertzian response on a specimen of elastic
modulus E = 400 kPa, while the cyan plot represents the force response on a
finite element (FE) model2 of a linearly elastic specimen of thickness h = 4µm
and same elastic modulus (E = 400 kPa). As expected, the force response
for the FE model is higher due to the underlying stiff substrate [44]. When
these force-indentation plots are recast into σ∗ - ∗ plots using Eqn. 5.3(a-c),
it becomes apparent that the FE solution shows significant deviation from the
2The details of the finite element (FE) model used to generate this plot are given in section 5.5.
113
Hertzian response [Fig. 5.4(b)], which is linear (from Eqn. 5.4). This nonlin-
earity occurs in the FE response not due to any hyperelasticity in the material,
but due to the incorrect formulation in the indentation strain ∗ in Eqn. 5.2(b)
which does not account for specimen depth. In fact, research groups have also
defined representative indentation strain as ∆/t and ∆/R [77]; however, there
is no general consensus on this matter. Apart from the aforementioned works,
0 200 400 600 800 1000
0
500
1000
1500
2000
2500
FE
Hertz
0 0.02 0.04 0.06 0.08 0.1 0.12
0
50
100
150
200
250
300
350
FE
Hertz
(a) (b)
Figure 5.4: Effect of finite deformation on the force response during spherical inden-
tation. (a) Simulated force response on a finite element model of a tissue specimen
(cyan) and a Hertzian force response of the same elastic modulus (purple) and (b)
corresponding indentation stress (σ∗) - indentation strain (∗) using Eqns. (5.3a -
c).
researchers have proposed inverse FE based methods to estimate hyperelas-
tic material parameters from AFM force curves [78][79][80]. Most of these
methods, however, do not consider the problem of finite-indentation.
• Finite Indentation & Nonlinear Elasticity (H3,H4,H5): Amongst the first
groups to study the coupled problem of finite indentation and material nonlin-
earity was that of Kosta and Yin [66]. Using finite element methods and the
114
“pointwise modulus” [1] approach, the authors studied the effects of indenter
geometry, hyperelasticity and sample thickness for AFM based characteriza-
tion studies. Likewise, Oommen and Val Vilet [77] studied the effects on
indentation on thin hyperelastic polymeric films using FE methods. However,
most of these approaches are primarily exploratory, and have little implications
in estimating the elastic modulus from any given AFM force-curve.
In this context, the work of Long et al. [59] gains prominence: Long et al. ex-
tended Dimitriadis’ approach [44] of appending correction factors to the Hertz
model, and both finite-indentation and material nonlinearity were accounted
for in [59]. Using a judicious choice of exponents of κ ( = R∆/h2, see Eqn.
5.2) in the correction terms, Long et al. fitted FE force-indentation data from
a neo-Hookean hyperelastic material to an analytical expression of force F as
a function of R,∆, κ (see Eqn. 4.7 for full expression), and obtained critical
bounds of R/h and ∆/h in which the expression was valid. Since the force
F could be explicitly related to ∆, it was also possible to use Long’s contact
model in the EIV-Bayesian Changepoint analysis for uncertainty quantifica-
tion in chapter 4.
However, as seen later in section 5.4, even neo-Hookean hyperelasticity is not
adequate to describe some of the force curves that are acquired during AFM
raster scanning on breast tissue specimens.
• Isotropic and Homogeneity (H6) Mechanical anisotropy has been shown to
exist in bone [81], skeletal muscle [82] and recently, vocal fold tissue [83]. In
115
[83], the authors carried out indentation studies on vocal fold tissue samples
sectioned at different depths along mutually perpendicular axes and the ratio
of their elastic modulus was used to quantity the degree of anisotropy.
Heterogeneity in the indentation direction has been addressed using numer-
ical methods. Unnikrishnan et al. [84] modeled inhomogeneity in cellular
cytoplasm by modeling the actin cortex as a hyperelastic material and inner
cytoplasm as a fiber-reinforced composite.
5.4 Problem Statement and Proposed Solution
From the preceding discussion in section 5.3 and the results in Fig. 5.3, it be-
comes clear that accounting for violations (H3-H6) in AFM data analysis is critically
important to ensure accurate mechanical characterization. In this context, Long’s
analytical relationship [59] appears particularly attractive since its accounts for vio-
lations (H3 - H5), however, it must be emphasized that neo-Hookean hyperelasticity
is a first-order approximation and it does not capture high degrees of nonlinearities
at large indentations. As a result, it might not be capable of capturing the spatially
heterogeneous tissue response as seen in Fig. 5.3.
This is well demonstrated in Figs. 5.5 and 5.6 where two experimental force
curves E1 and E2 are analyzed by implementing Hertz, Dimitriadis’s and Long’s
contact models in the EIV-Bayesian Changepoint algorithm described in chapter
4. In force curve E1 (Fig. 5.5), it is evident that all three contact models fit the
116
Table 5.1: Summary of implementation of the EIV-Bayesian Changepoint algorithm on two experi-
mental force curves E1 and E2 (Figs. 5.5 and 5.6 respectively acquired on a tissue specimen.)
Elastic Modulus (kPa) Contact Point Est. Indentation R2
Mean 95% C.I. Mean 95% C.I. ∆ˆ = δn − δbkˆc
Force Hertz 178.30 [153.19 205.98] 447 [446 449] 461.5 nm 0.999
Curve (E1) Dimitriadis 110.48 [95.15 122.98] 427 [425 428] 490.7 nm 0.999
(Fig. 5.5) Long 119.01 [105.41 131.98] 433 [429 434] 481.8 nm 0.999
Force Hertz 151.83 [132.31 175.65] 141 [139, 143] 433.1 nm 0.908
Curve (E2) Dimitriadis 106.27 [92.78 122.92] 141 [139, 143] 433.1 nm 0.926
(Fig. 5.6) Long 114.382 [99.76, 132.27] 141 [139, 143] 433.1 nm 0.917
117
deflection data well (overall R2 > 0.999)3. The Hertzian elastic modulus is under-
standably higher compared to the other two due to infinitesimal assumptions. The
estimated elastic modulus using both Dimitriadis’s and Long’s model are close to
each other, which can be explained by observing that the Long’s model reproduces
the result from Dimitriadis’s model for indentation regimes ∆/h close to 0.1 [59]
(The neo-Hookean model is “nearly” elastic for low deformation regimes).
Fig. 5.6 shows the implementation on force curve E2 acquired in the spatial neigh-
0 200 400 600 800
−10
0
10
20
30
40
50
0 200 400 600 800
−10
0
10
20
30
40
50
0 200 400 600 800
−10
0
10
20
30
40
50
420 440 460 480 500
−2
0
2
4
6
420 440 460 480 500
−2
0
2
4
6
420 440 460 480 500
−2
0
2
4
6
(a) (c)(b)
Figure 5.5: Results of Gibbs Sampling on the posterior distribution using (a) Hertz
contact model, (b) Dimitriadis’s contact model and (c) Long’s contact model on
force curve E1.
borhood of force curve E1. It can be observed that none of the aforementioned
contact models produce a good fit (overall R2 ∼ 0.90− 0.92). Typically, such force
curves are abundantly acquired during raster scans on tissue specimens, where a
single contact model does not appropriately fit all force-curves acquired on a tissue
region-of-interest. The poor fit quality indicates that even a neo-Hookean hypere-
lastic constitutive model is not sufficiently nonlinear enough to explain the sharp
tissue response in force curves such as those shown in Fig. 5.6. Furthermore, it
3By the term “overall R2”, the fit quality for both non-contact and contact regimes are meant.
118
0 100 200 300 400 500
−10
0
10
20
30
40
50
60
0 100 200 300 400 500
−10
0
10
20
30
40
50
60
0 100 200 300 400 500
−10
0
10
20
30
40
50
60
(a) (c)(b)
Figure 5.6: Results of Gibbs Sampling on the posterior distribution using (a) Hertz
contact model, (b) Dimitriadis’s contact model and (c) Long’s contact model on
force curve E2.
is worthwhile to note that both force curves registered similar indentation range
(∆ˆ/h ∼ 0.115 − 0.122 for force curve E1 and ∆ˆ/h ∼ 0.108 for curve E2 for a
tissue specimen 4 µm thick) though their fit quality were vastly different. This in-
dicates that large indentation is not the only cause of the observed tissue response
in Fig. 5.6, were it to be so, a 10% sample indentation in force curve E1 too would
have shown a poor fit quality. This is indicative of the fact that the tissue specimen
was inherently nonlinear beyond the neo-Hookean approximation in force curve E2.
The immediate implication of this observation is that truncation of force-indentation
data, as is often used to limit the fitting range to the linear elasticity regime [78],
might still lead to inaccurate estimation of the elastic modulus if the underlying
constitutive material model for the force curve is not established.
Up to this point, the effect of the contact point has not been discussed in
estimating hyperelastic constitutive parameters from AFM force curves. In chapter
4, the EIV-Bayesian Changepoint algorithm that was used for contact estimation
had one critical assumption: the contact model for an AFM force curve was known
119
beforehand. From the preceding discussion, it is clear that Long’s contact model
(with neo-Hookean hyperelasticiy) would not be able describe force curves such as
the one shown in Fig. 5.6. Therefore, the contact estimation procedure would need
to be re-validated for force curves that show high degree of material nonlinearity.
Hypothesis: By incorporating finite element modeling of spherical inden-
tation of nonlinear hyperelastic tissue response in the mechanical characterization
process, it is possible to improve accuracy and robustness in AFM indentation data
analysis.
To address the problem of material nonlinearity beyond neo-Hookean approx-
imation and its associated contact estimation problem, the following solution is
proposed:
• FE analysis of hyperelastic tissue specimens under AFM loading (Assumptions
H3-H5:) To our knowledge, there are no existing explicit contact models for
capturing finite-indentation and material nonlinearity simultaneously beyond
neo-Hookean hyperelasticity. This necessitates use of numerical techniques to
accommodate data analysis for highly nonlinear force curves and to this end,
2D axisymmetric finite element models of tissue specimens with an exponential
form of the strain energy density function have been developed that sufficiently
captures the nonlinear response as evidenced by Fig. 5.6. The details of the
FE modeling is described in section 5.5.
• Contact Estimation: A weighted least-squares based approach is proposed to
estimate the contact point from force curves with high degrees of nonlinearity.
120
Details of the contact estimation algorithm are discussed in section 5.6.
It should be mentioned here that during AFM raster scanning on breast tis-
sue specimens, force curves such as Fig. 5.7 were frequently encountered that indi-
cated that the assumption of homogeneity in the indentation direction was violated
[Hertzian assumption (H6)]. These force curves were excluded from the data analysis
outlined in this chapter. Moreover, isotropic assumptions in the tissue for modeling
and analysis were also retained in this chapter.
−0.5 0 0.5
0
50
100
150
Z−position, z, (µm)
De
fle
ct
io
n,
 d
 (n
m
)
0 100 200 300 400 500
0
20
40
60
80
Indentation, ∆ (nm)
Fo
rc
e,
 F
 (n
N)
(b)(a)
Figure 5.7: Evidence of lack of material homogeneity in the direction of indentation.
(a) AFM force curves. Please note that the force curves have been displaced to have
the same approximate contact point at (0,0).(b) Corresponding force-indentation
curves. Note the change in the nature of the force-indentation curves at 200 nm and
300 nm indentation.
5.5 Finite Element Modeling
5.5.1 Discretization and Parameter Estimation
The FE analysis was carried out in ABAQUS v6.8 (Dassault Systemes, RI,
USA). Taking advantage of the symmetry of the indentation process, a 2D axisym-
121
metric model of the indenter and the tissue specimen was used. The quarter of the
spherical indenter in contact with the specimen was modeled as an axisymmetric
rigid indenter.
The tissue specimen consists of 135 axisymmetric, quadrilateral hybrid ele-
ments (CAX4H) (15 radially, 9 axially). The length of the tissue model was 6 µm.
A fine mesh was used near the indentation region, with a coarser mesh away from it
(see Fig. 5.8). No friction was assumed between the indenter and the tissue specimen
model. The nodes at the interface of the specimen and the indenter were considered
traction free, while the nodal displacements in the radial and axial directions were
constrained. Nodes on the axis were constrained to move only along the direction of
indentation. Quasistatic equilibrium solutions were acquired for incremental axial
displacements of the indenter. In a mesh convergence study, it was verified that
mesh size or the length of the tissue model did not have significant impact on the
analysis results.
With any single or multi-parameter constitutive model whose explicit force-
indentation relationship is unknown, estimating the material properties typically
involves solution of an inverse-finite element problem where the simulated response
is curve-fitted to experimental data. This is obtained by minimizing the squared
norm of the error and is given by the following expression:
argmin
θ
∑
||Fexp − FFE(θ)||2 (5.5)
where Fexp and FFE represent the experimental and simulated force response re-
122
(a)
(b)
Figure 5.8: Finite element models: (a) undeformed mesh, (b) deformed mesh.
123
spectively, and θ represents a vector of material parameters defined in the consti-
tutive material model. The minimization problem was implemented in MATLAB
(Mathworks, Inc) using the fminsearch routine which is a heuristic optimizer using
the Nelder-Mead simplex algorithm [85]. The optimizer iterates from an arbitrary
initial point towards the optimum solution until the incremental reduction in the
objective function is lesser than a predefined tolerance. A flowchart indicating the
steps of the inverse FE approach is shown in Fig. 5.9.
Finite Element Analysis
Update 
Initialize 
Equalize
and       
    
Has the solution 
converged ?
Obtain 
End
 
Yes
No
Figure 5.9: Flowchart of the algorithm used to obtain optimized material properties
using the inverse FE approach.
124
Note that in Fig. 5.9, the block “contact estimation” has been included in the
flowchart without explicitly describing a technique for it. The contact estimation
procedure is discussed in section 5.6.
5.5.2 Constitutive Material Models
The tissue specimen is modeled as an isotropic, incompressible and hypere-
lastic material. Isotropic assumptions in the indentation direction were considered
valid because datasets such as those in Fig. 5.7 were not included in this study.
Additionally incompressibility conditions were assumed in the tissue specimens due
to the hydrated nature of the tissue specimens in phosphate buffered saline (PBS)
during AFM indentation experiments.
In order to model force curves like force curve E2 shown in Fig. 5.6, significant
material nonlinearity needs to be considered in the constitutive material model. The
force response for the following two incompressible hyperelastic materials character-
ized by the following strain energy functions were examined.
WY EOH = C1(I1 − 3) + C2(I1 − 3)2 + C3(I1 − 3)3 (5.6a)
WEXP =
C
b [e
b(I1−3) − 1] (5.6b)
where (C1, C2, C3) and (C, b) are hyperelastic material parameters of the Yeoh and
Fung hyperelastic models respectively and I1 is the first invariant of the Cauchy-
Green deformation tensor, expressed in terms of the principal stretch ratio λ1, λ2, λ3
125
as:
I1 = λ21 + λ22 + λ23 (5.7)
where λ1λ2λ3 = 1Eqn. (5.6a) represents the Yeoh hyperelastic material [86] which
can be viewed as a polynomial extension of the neo-Hookean model (with C2, C3 =
0), and is commonly used to model rubber elasticity. Eqn. (5.6b) represents an
exponential form proposed by Fung [87], which has been previously used to model
biomaterials. The presence of higher orders coefficients such as (C2, C3) and b in the
Yeoh and Fung models respectively makes these two models particularly attractive in
describing highly nonlinear force curves such as Fig. 5.6. While the Yeoh material
model was pre-defined in the ABAQUS material library, Fung’s model was not.
Consequently, the user subroutine UHYPER was employed to specify the strain
energy function denoted by Eqn. (5.6b) for the FE analysis.
The hyperelastic coefficients can be related to the specimen’s materials prop-
erties using the relation:
µ0,Y EOH = 2
∂WY EOH
∂I1
∣∣∣∣
I1=3
= 2C1 (5.8a)
µ0,EXP = 2
∂WEXP
∂I1
∣∣∣∣
I1=3
= 2C (5.8b)
where µ0 is the initial shear modulus, which is related to the initial elastic modulus
E0 by:
E0 = 2(1 + ν)µ0 (5.9)
When incompressibility conditions (ν = 0.499) apply, E0 = 6C1 and E0 = 6C for
126
Yeoh and Fung material models respectively.
0 100 200 300 400 500
−50
0
50
100
150
200
0 100 200 300 400 500
−50
0
50
100
150
200
(a) (b)
Figure 5.10: Performance of the (a) Yeoh and (b) Fung hyperelastic models in
describing the force response for force curve shown in Fig. 5.6.
Fig. 5.10 shows the force curve E2 from Fig. 5.6 overlaid with FE simula-
tion results of indentation on a Yeoh [Fig. 5.10(a)] and Fung hyperelastic material
[Fig. 5.10(b)]. The contact point was assumed to be the bkˆc = 141, as obtained
from the EIV-Bayesian Changepoint algorithm, leading to an indentation of 433.1
nm (from Table 5.1) in the simulation. The results demonstrate the applicability of
higher order hyperelastic material models to reproduce force curves with pronounced
nonlinearity (overall R2 > 0.99). The estimated hyperelastic parameters using an
inverse FE analysis were (C1, C2, C3) = (5.94, 134.72, 4900) kPa and (C, b) = (7.55
kPa,45.0) for the Yeoh and Fung hyperelastic models respectively. Using Eqn. (5.8),
the initial elastic modulus was computed to be E0 = 35.68 for Yeoh’s model and
E0 = 45.33 kPa for Fung’s model.
Several conclusions can be drawn from the results in Fig. 5.10. Firstly, esti-
mates of higher order coefficients in the Yeoh model (C2, C3) are significantly greater
127
than C1, thereby demonstrating that sharp material nonlinearity in the force curve
can be successfully captured by (C2, C3), while C1 quantifies only the initial tissue re-
sponse to indentation. Compared to the case where only a single material parameter
is extracted from the force curve, i.e. the case when C2 = 0, C3 = 0: a neo-Hookean
hyperelastic model where force and indentation explicitly relates by Long’s contact
model [59], the elastic modulus is significantly overestimated (114.382 kPa, see ta-
ble 5.1), the obvious implication being that more than one material parameter is
necessary to describe the force curves such as the one in Fig. 5.10 completely.
Secondly, the quality of fit and the estimated initial elastic modulus being
similar for both the Yeoh and Fung hyperelastic models, a pertinent question comes
up: which hyperelastic model should be selected for modeling force curves shown in
Fig. 5.10 ?
A parametric study is demonstrated in Fig. 5.11, where the hyperelastic pa-
rameters for the Yeoh and Fung models ([Fig. 5.11(a-c)] and [Fig. 5.11(d,e)] respec-
tively) are varied individually, and their corresponding force-indentation curves are
obtained. Both the Yeoh and Fung models utilize a single parameter, C1 and C,
respectively to capture the initial tissue response. For increasing material nonlinear-
ity, Yeoh’s hyperelastic model uses two parameters, C1 and C2, while Fung’s model
uses only one, b. In this context, the Fung model is appears particularly attractive
because of one less hyperelastic parameter compared to the Yeoh model. From the
perspective of solving an inverse-FE problem (typically for estimating material prop-
erties from an experimental AFM force curve), it is also likely that the Fung model
would pose less existence and uniqueness problems compared to the Yeoh model,
128
0 100 200 300 400 500
0
50
100
150
200
0 100 200 300 400 500
0
50
100
150
200
0 100 200 300 400 500
0
100
200
300
400
500
0 100 200 300 400 500
0
100
200
300
400
500
0 100 200 300 400 500
0
100
200
300
400
500
(a) (b) (c)
(d) (e)
Figure 5.11: Effect of varying the hyperelastic parameters of the Yeoh (a,b,c) and
Fung (d,e) hyperelastic models individually.
where material nonlinearity is coupled in between C2 and C3. In contrast, Fung’s
model is fairly decoupled: C quantifies the initial elastic modulus, while b describes
the degree of material nonlinearity. For these reasons, the Fung hyperelastic model
is used to describe the observed hyperelasticity in the AFM force curves.
5.5.3 Sensitivity Analysis
Since the Fung model captures material nonlinearity in AFM force curves
through the parameter b, it is instructive to examine the sensitivity of b to varying
degrees of nonlinear elasticity in AFM force curves. The estimated value of b (=
45.0) in Fig. 5.10 raises an interesting question: is there an existence of a lower
bound of b which corresponds to a linear elastic response (the analytical analogue of
129
which is Dimitriadis’s [44] model) ? If this is the case, then the parameter b could
be used to potentially differentiate between linear and nonlinear elastic responses
in the tissue specimen, which could serve as an additional mechanical signature for
tissue health.
5
10
15
20
25
30
300350400450500
5
10
15
20
25
30
300350400450500
5
10
15
20
25
30
300350400450500 300
0.995
0.996
0.997
0.998
0.999
1
5
10
15
20
25
30
300350400450500
5
10
15
20
25
30
300350400450500
5
10
15
20
25
30
300350400450500 300
0.995
0.996
0.997
0.998
0.999
1
0
1
2
3
4
5
6
300350400450500
0
1
2
3
4
5
6
300350400450500
0
1
2
3
4
5
6
300350400450500 300
0.995
0.996
0.997
0.998
0.999
1
(a) (b) (c)
(d) (e) (f)
(g) (h) (i)
Figure 5.12: Sensitivity of the Fung model parameters at varying indentations.
(a)∆ = 100 nm, (b) ∆ = 200 nm, (c) ∆ = 400 nm, (d) ∆ = 600 nm, (e) ∆ = 900
nm and (f) ∆ = 1400 nm. Also shown in the third row (g-i) are enlarged figures of
(d-f) where the range of b is reduced to 1 − 6. The white dotted line at 400 kPa
indicates the elastic modulus of the linearly elastic material.
A sensitivity analysis was performed to examine values of E0(= 6C) and b
that resulted in a good fit with a F − ∆ curve obtained from a linearly elastic
material, and the results are shown in Fig. 5.12. The target F − ∆ curve was
130
obtained by varying the indentation from ∆ = 100 nm (∆/h = 0.025) to ∆ = 1400
nm (∆/h = 0.35) on a linearly elastic specimen of elastic modulus 400 kPa and 4
µm thickness.
The ranges of E0 and b that fit well with the linearly elastic force response
can be compactly expressed by the region S(∆,R2FE,thresh) ∈ R
2 where S(∆,R2FE,thresh) =
{(E0, b)|R2∆(E0, b) > R2FE,thresh}, where R2FE,thresh indicates a pre-specified threshold
R2 value. The term R2FE,thresh has implications for an inverse FE analysis during
parameter estimation - it determines the termination criteria for an inverse FE
solution. In Fig. 5.12, R2FE,thresh was chosen to be 0.995 (see colorbar in Fig. 5.12)4.
For low indentation regimes such as ∆ = 100 nm, [Fig. 5.12(a)], moderately
hyperelastic specimens (b ∼ 25 − 30)5 produce a linearly elastic response. With
increasing ∆, the domain of S(∆,R2FE,thresh) progressively shrinks. The immediate
implication of this observation is that parameter estimation using an inverse FE
approach is likely to yield inaccurate estimates of b for shallow indentation regimes
due to the relatively large domain of S(∆,R2FE,thresh), and it is unlikely that it would
be possible to differentiate between linear and hyperelastic tissue response for these
regimes. This effect can be avoided during AFM experiments by using a relatively
stiff cantilever and specifying (1) a moderately large force setpoint in force-controlled
mode, or (2) by setting a large indentation setpoint in indentation-controlled mode.
However, for larger indentations, the estimated values of b can be used for
differentiating between linear and hyperelastic tissue response. This is evident by
4Minimizing the squared norm of the residuals is equivalent to maximizing the R2.
5The peak b at ∆ = 100nm was 67.0 for R2FE,thresh = 0.995 (not shown).
131
observing that for an indentation of ∆ = 433 nm in force curve E2 of Fig. 5.10 (see
Table 5.1), the estimated b was 45.0, whereas the b is upper-bounded in S(400,0.995)
by a value of 10-10.5 [see Fig. 5.12(c)]. It now becomes possible to construct a clas-
sification curve (see Fig. 5.13), where the linearly elastic tissue response is bounded
by bmaxR2FE,thresh(∆), which is defined as:
bmaxR2FE,thresh(∆) := max{b|b ∈ S(∆,R2FE,thresh)} (5.10)
0 200 400 600 800 1000 1200 1400
0
20
40
60
80
100
120
Estimated 
from Fig. 5.10
Figure 5.13: Classification curve, showing the variation of bmaxR2FE,thresh(∆) with ∆ for
different values of R2FE,thresh. Also shown are the curve fits for bmaxR2FE,thresh(∆).
132
5.6 Contact Estimation
In chapter 4, the problem of contact estimation associated with AFM force
curves on soft tissue was discussed and a probabilistic EIV-Bayesian Changepoint
algorithm was proposed for uncertainty quantification in AFMmechanical character-
ization. One critical assumption was made in EIV-Bayesian Changepoint algorithm:
the post-contact F −∆ analytical model was assumed known. It might be recalled
that the EIV-Bayesian Changepoint analysis was restricted to those force curves
only for which the post-contact regime could be described by the analytical contact
model used. For force curves with pronounced hyperelasticity where no contact
model exists, the contact estimation problem needs to be revisited.
It was initially envisaged that contact estimation could be integrated in the
inverse-finite element algorithm. This would obviate the need to have a explicit
analytical contact model for estimating the contact point. Essentially, this approach
involved augmenting the cost function (from Eqn. 5.5) with the unknown contact
point index k and simultaneously estimating the set of parameters (E0, b, k) that
minimized the augmented cost, as shown below:
arg min
E0,b,k
k∑
i=1
[Fi − β11 − β12δi]2 +
n∑
i=k+1
[Fi − FFE(E0, b)]2 (5.11)
where β11, β12 are pre-contact regression parameters, which could be estimated by
traditional least squares fitting. It might be recalled from chapter 4 that δ, was
obtained by subtracting the AFM probe, d, from the vertical probe position, z, i.e.,
133
δ = z − d.
However, several numerical difficulties were encountered with this approach.
Standard gradient-based algorithms and heuristic approaches such as the Nelder-
Mead simplex algorithm [85] converged to local minimas, primarily because the
curvature of the cost function was significantly lower along the direction of contact
(i.e k) for smooth transitions compared to the directions of E0 and b. Since an
inverse-FE approach by itself is computationally expensive, which might have to
be implemented on a large number force curves on a given tissue region of interest
(ROI), the approach of simultaneous estimation of the contact point and material
parameters was not pursued any further.
In the subsequent section, two approaches for estimating the contact are de-
scribed. First, the bidomain polynomial (BDP) approach without continuity con-
straint is discussed, which was proposed by [1]. This can be viewed as a determin-
istic version of a Bayesian Changepoint problem6 that was discussed in chapter 4,
nonetheless, a description of this approach is presented for the sake of completeness.
Then, a continuity enforced weighted BDP approach is discussed. It should be noted
that both these approaches require an explicit contact model.
5.6.1 BDP approach without continuity [1]
It might be recalled from chapter 4 that the contact estimation in an AFM
force curve was cast as a two-regime regression model with a unknown changepoint
6Since spring constant variations are not considered in this section, the EIV-Bayesian Change-
point algorithm reduces to a regular Bayesian Changepoint algorithm.
134
separating the two regimes. The two-regime regression model (Eqn. 4.12 from
chapter 4) is restated below:
Fi =



β11 + δiβ12 + 1 if i ≤ k
β21 + fiβ22 + 2 if k + 1 ≤ i ≤ n
(5.12)
where Fi is the observed force obtained by multiplying the deflection di by the
probe spring constant kc, β1 = [β11 β12]T and β2 = [β21 β22]T are the pre-
contact and post-contact regression coefficients respectively, δi is the transformed
variable obtained as δi = zi − di, k ∈ (1, n) is the unknown contact index for n
datapoints, 1 iid∼ N(0, σ12), 2 iid∼ N(0, σ22) are the error terms in the pre-contact
and post-contact regimes respectively, and fi is the contact model.
The bidomain polynomial (BDP) approach seeks to obtain optimal estimates
(kˆ, βˆ1, βˆ2) that minimizes the error between the observed force readings Fi and the
fitted force estimates from the regression model of Eqn (5.12).
As in chapter 4, the force data can be written in compact form as F k1 =
[
F1 F2 . . . Fk
]T
, F k2 =
[
Fk+1 Fk+2 . . . Fn
]T
and F =
[
F k1 F k2
]
. The
right hand side of Eqn. (5.11) form the design matrices given by:
Xk1 =


1 1 . . . 1
δ1 δ2 . . . δk


T
135
Xk2 =


1 1 . . . 1
fk+1 fk+2 . . . fn


T
(5.13)
The optimal estimate for the regression coefficients (βˆk1, βˆ
k
2) at each k is there-
fore given by:
(βˆk1, βˆ
k
2) = arg minβ1,β2
||F k1 −Xk1β1||2 + ||F k2 −Xk2β2||2 (5.14)
Eqn. (5.14) is the least squares problem, whose solution is given by:
βˆk1 = [(Xk1)TXk1]−1(Xk1)TF 1
βˆk2 = [(Xk2)TXk2]−1(Xk2)TF 2 (5.15)
The BDP estimate of the contact point obtained by performing a line search on the
sum of squared residuals (SSR), and is given by:
kˆ = arg min
1<k<n
SSR (5.16)
where
SSR = ||F k1 −Xk1βˆ
k
1||2 + ||F k2 −Xk2βˆ
k
2||2 (5.17)
5.6.2 Continuity-constrained weighted BDP approach
In the unconstrained BDP approach discussed in section 5.6.1[1], no continuity
is considered at the changepoint, i.e. β21 is not necessarily equal to β11 + δkβ12.
136
Not enforcing continuity did not significantly impact the contact point estimate in
chapter 47, primarily because the contact model fitted the post-contact deflection
data well. For sharply nonlinear force curves for which no contact model exists, it
is unlikely that the contact estimation would be accurate.
As a result, a modification to the BDP approach is proposed, with continuity
enforced at the contact point. This leads to the following constraint:
β21 = β11 + δkβ12 (5.18)
Due to the constraint in Eqn. (5.18), the regression parameter space reduces to
three unknowns: β˜ = [β11 β12 β22]T . The modified regression model from Eqn.
(5.12) are given by the following equations:


F1
...
Fk
Fk+1
...
Fn


=


1 δ1 0
... ... ...
1 δk 0
1 δk fk+1
... ... ...
1 δk fn


×


β11
β12
β22


(5.19)
7It might be recalled that the EIV-Bayesian Changepoint implementation on simulated data
(section 4.4.2) yielded the posterior mean of the contact point within 3 datapoints of the ground
truth.
137
In vector form, these equations are given as:
F︸︷︷︸
n × 1
= X˜k︸︷︷︸
n × 3
β˜k︸︷︷︸
3 × 1
Since the contact model fi is unlikely to fit the entire post-contact force data for
highly nonlinear force curves, weights are added to the residuals of the regression
model to restrict the fitting range near the contact transition. This is accomplished
by defining a diagonal weight matrix W k = diag(
√
wk11,
√
wk22, . . . ,
√
wknn) where
W k ∈ Rn×n.
The weights are assumed Gaussian, and are given as follows:
wii = w0e
− 12σ2w
(i−k)2 (5.20)
Note the dependence of the weights around the candidate contact point k -
this ensures that the contact transition is captured by the regression model.
Similar to the unconstrained BDP approach, the optimal regression coefficient
estimate ˆ˜βk for a given contact point k is given by:
ˆ˜βk = argmin
β˜
||W k(F − X˜kβ˜)||2 (5.21)
The weighted least squares solution ˆ˜βk is given by:
ˆ˜βk = [(X˜k)T (W k)TW kX˜k]−1(X˜k)TW kF (5.22)
138
The constrained BDP estimate of the contact point obtained by performing a
line search on the sum of squared residuals (SSR) and is given by:
kˆ = arg min
1<k<n
SSR (5.23)
where
SSR = ||(F − X˜k ˆ˜βk)||2 (5.24)
5.6.3 Choice of fi for automated data analysis
Central to both BDP model (with and without continuity constraint) is the
post-contact regression model denoted by fi. As discussed previously, analytical
force-indentation relationships are not available beyond neo-Hookean hyperelastic-
ity. In this context, it must be noted that the need to estimate the contact point
and estimate hyperelastic material parameters (using inverse FE) arises only when
there is prominent hyperelasticity in the force curves - otherwise, the proposed ap-
proach in chapter 4 suffices in accurately quantifying mechanical properties from
AFM force curves.
It therefore seems prudent to use Dimitriadis’s force-indentation relationship
as fi in the BDP approaches suggested in sections 5.6.1 and 5.6.2. Using Dimitri-
139
adis’s model [44], fi is given by:
fi =



0 if i ≤ k
(δi − δk)1.5 + 1.133
√γ(δi − δk)2 + 1.283γ(δi − δk)2.5 if k + 1 ≤ i ≤ n
+0.769γ√γ(δi − δk)3 + 0.0975γ2(δi − δk)3.5
(5.25)
where γ = R/h2. If Dimitriadis’s model fit (R2Dim) exceeds a threshold (R2Dim,thresh)
after the contact point is established either using the regular BDP or the con-
strained weighted BDP approach, no inverse FE analysis is required since the force
curve represents linear elasticity in the material, which is adequately captured using
Dimitriadis’s model. In the case where the post-contact fit to Dimitriadis’s model is
poor, inverse FE estimation is invoked with the Dimitriadis’s solution as the initial
condition, which is terminated when the R2FE reaches the threshold R2FE,thresh. A
pseudocode describing the proposed automation strategy is given in below (Algo-
rithm 1).
5.7 Results
In this section, the algorithm shown in Algorithm 1 is implemented on AFM
force curves of varying degrees of elasticity. The proposed approach is first validated
on simulated force curves. Then, the performance of the proposed approach is
demonstrated on a tissue elastic map.
140
input : AFM force curve, R2Dim,thesh,R2FE,thresh.
output: Material parameters.
Multiply deflection with probe spring constant to obtain force;
Contact Estimation: Estimate contact point using BDP based
approaches;
Dimitriadis Curve fitting: Fit F −∆ to Dimitriadis’s contact model
to obtain EDim, R2Dim;
if R2Dim > R2Dim,thesh then
Linear Elasticity: Estimated material parameter: Eˆ ← EDim;
else
Hyperelasticity:
end
initialize C0 ← EDim/6 , b← 30;
Invoke Inverse FE routine with R2FE,thresh;
Estimated material parameters: Eˆ0 ← 6Cˆ0 , bˆ;
Algorithm 1: Pseudocode for automated force curve analyis
5.7.1 Parameter reidentification on simulated force curves
Three AFM force curves (S1-S3) are simulated with varying degrees of hyper-
elasticity. The parameters used to construct force curves are given in Table 5.2.
Additionally, the force values were perturbed by Gaussian noise [∼ N(0, 0.05)]. The
threshold parameters R2Dim,thresh and R2FE,thresh were set to 0.980 and 0.999 respec-
tively.
Table 5.2: Parameters used to construct the simulated force curves.
(Contact k, Inden- Linear Hyperelastic
Datapoints n) tation Elastic Modulus Parameters (E0, b)
Force Curve S1 (501,900) 400 nm 200 kPa ×
Force Curve S2 (501,900) 400 nm × (65 kPa,60)
Force Curve S2 (501,900) 400 nm × (45 kPa,120)
Performance of the unconstrained BDP [1] and continuity-constrained weighted
BDP methods applied to force curves S1, S2 and S3 are shown in Fig. 5.14. The
141
weighting parameters used were w0 = 30, σw = 3. In addition, the performance of
the continuity-constrained weighted BDP is investigated without additional weights,
i.e. W k = I, where I is the identity matrix. It is clear from Fig. 5.14(a-c) that
when the material is linearly elastic, all three approaches produce close contact es-
timates compared to the ground truth (501th datapoint). This is primarily because
the Dimitriadis contact model, i.e. [44] produces an elastic response - which was
used to obtain force curve S1. The unconstrained BDP estimate [Fig.5.14(a)] can
be further corrected using a Bayesian Changepoint scheme, (as demonstrated in
chapter 4) and it is likely that the 95% confidence intervals will contain the ground
truth.
With increasing hyperelasticity in force curves S2 and S3, it is evident that the
performance of all three approaches degrade. However the continuity-constrained
weighted BDP approach significantly outperforms the other two approaches. This
is primarily because the fitting range is restricted by the assigned weights. In the
absence of any weights, the continuity constrained BDP still outperforms the un-
constrained BDP approach, however the estimation error is significant.
One of the shortcomings of the continuity-constrained weighted BDP approach
is that weights require to be tuned to produce correct estimates of the contact point
for sharply nonlinear force curves. It is possible to reduce the contact estimation
error in Fig. 5.14(h) further by increasing the weights, however, very large weights
(σi < 3) tends to yield false minimas in the sum of squared residuals (SSR) function.
Typically, w0 = 30, σw = 3 works optimally for force curves with moderate to large
nonlinear elasticity (such as force curve S2) in the force curves. For force curves
142
such as S3 with very large b, however, this contact estimation error in inevitable
using the current approaches.
 S2
Ground truth
 S3
 S1
Contact estimate with unconstrained BDP 
200 300 400 500 600 700 800
0
2
4
6
8
x 107
0
2
4
6
8
x 105
200 300 400 500 600 700 800
0
0.5
1
1.5
2
x 105
300 400 500 600 700 800 900
0
x 107
0
10
x 105
300 400 500 600 700 800 900
0
5
10
15
x 105
300 400 500 600 700 800 900
0
x 107
0
2
4
6
8
10
x 106
300 400 500 600 700 800 900
0
2
4
6
8
10
12
x 106
0 200 400 600 800
−50
0
50
100
150
200
0 200 400 600 800
−50
0
50
100
150
200
250
300
0 200 400 600 800
0
200
400
600
800
Datapoints
Datapoints
Datapoints Datapoints
Datapoints
DatapointsDatapointsDatapoints
Datapoints
SS
R
SS
R
SS
R
SS
R
SS
R
SS
R
SS
R
SS
R
SS
R
Fo
rc
e 
(n
N)
Fo
rc
e 
(n
N)
Fo
rc
e 
(n
N)
Sum of Squared Residuals (SSR) 
using unconstrationed BDP
Contact estimate with 
weighted constrained BDP
Sum of Squared Residuals (SSR) 
using weighted constrained BDP
Contact estimate with unweighted 
 constrained BDP 
Sum of Squared Residuals (SSR) 
using unweighted constrained BDP
(a) (b) (c)
(f)(e)(d)
(g) (h) (i)
516 501 501
687 659
517
593 766780
Figure 5.14: Performance of the three contact estimation methods for simulated
force curves S1 - S3 (shown row-wise). The first column shows the sum of squared
residuals (SSR) using unconstrained BDP (red) for force curves S1-S3 (a),(d),(g) re-
spectively. The second column shows SSR variations using constrained BDP (blue
and green) for S1-S3 (b),(e),(h) respectively. The third column shows the corre-
sponding contact point on the force curves S1-S3 (c),(f),(i) respectively. The contact
point estimates using the unconstrained and constrained BDP approaches are also
shown at the base of the SSR curve.
Fig. 5.15 shows the Dimitriadis curve fitting and inverse FE solution steps
of the proposed automation strategy. The continuity-constrained weighted BDP
method is chosen for contact estimation since it outperforms the unconstrained BDP
143
and the unweighted continuity-constrained BDP. For force curve S1, the Dimitriadis
contact model produces a very good fit and therefore invocation of the inverse FE
routine is not necessary. For force curves S2 and S3, the inverse FE routine is
invoked. The estimated E0 and b from the inverse FE are very close to the ground
truth values for force curve S2. For S3 however, the there is an estimation error of
100%, which is a direct consequence of the contact estimation error.
0 100 200 300 400
0
50
100
150
200
0 100 200 300 400
0
50
100
150
200
250
300
0 100 200 300 400
0
50
100
150
200
250
300
0 50 100 150 200 250 300 350
0
200
400
600
800
0 50 100 150 200 250 300 350
0
200
400
600
800
S2
S3
S1
0 200 400 600 800
−50
0
50
100
150
200
0 200 400 600 800
−50
0
50
100
150
200
250
300
0 200 400 600 800
0
200
400
600
800
(a) (b)
(c) (d) (e)
(g)
501
517
593
501
501
501
(f) (h)
Figure 5.15: Implementation of the proposed algorithm for automated characteriza-
tion for simulated force curves S1 - S3 (shown row-wise). Note that the inverse FE
correction step was not required for curve S1.
144
The estimated parameters for all three force curves S1-S3 are tabulated in
Table 5.3.
Table 5.3: Summary of algorithm performance on force curves S1-S3. Units of E and
E0 are in kPa. G.T. indicates ground truth.
Force Contact Point k Material Parameters
Curve G.T. (k) (kˆ) G.T. (E) (E¯Dim) G.T. (E0) (Eˆ0) G.T. (b) (bˆ)
S1 501 501 200 196.97 × × × ×
S2 501 517 × × 65 71.82 60 60.78
S3 501 593 × × 45 94.54 120 150.69
5.7.2 Implementation on elastic map
In Fig. 5.3, an elastic map on a tissue ROI was shown where the spatial distri-
bution of the fit quality (R2) was not uniform throughout. The proposed algorithm
for automated characterization (Algorithm 1) is implemented on all the force curves
acquired in the elastic map to investigate if the inverse-FE correction yields more
accurate characterization results. Due to the computational cost involved for the
inverse-FE routine, a modest termination criteria for the inverse-FE procedure was
set (R2FE,thresh = 0.992), while R2Dim,thresh = 0.980 was used for detecting linear
elasticity in the force curve.
In Fig. 5.16, the implementation of the proposed algorithm is shown. Fig. 5.16(a)
shows the elastic modulus when no inverse FE corrections are implemented, i.e. the
Dimitriadis model is used to fit all the force curves after the contact point is esti-
mated using the continuity-constrained weighted BDP approach. The spatial dis-
tribution of the corresponding (R2) in Fig. 5.16(b) shows that the fit quality is still
145
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
0
100
200
300
400
500
600
700
800
900
1000 kPa
0
100
200
300
400
500
600
700
800
900
1000 kPa
(a) (b)
(c) (d)
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Figure 5.16: Implementation of the proposed algorithm for automated characteriza-
tion on a tissue ROI shown previously in Fig. 5.3. (a) Elastic map using the Dimitri-
adis model only without any Inverse FE corrections, (b) Corresponding post-contact
R2Dim fit, (c) Elastic map (Eˆ) with inverse-FE correction, (d) Corresponding R2net
fit, where R2net = max(R2Dim, R2FE). White spaces indicate force curves that were
not analyzed due to presence of multiple potential contact points such as those in
Fig. 5.7.
146
non-uniform. When the inverse FE correction is implemented, it becomes apparent
that the (R2) distribution improves significantly, because force curves with material
nonlinearity are now accounted for using Fung’s hyperelastic parameters.
In the elastic map, there were a few force curves that displayed multiple contact
points, as shown previously in Fig. 5.7. These force curves were not analyzed and
their corresponding pixel location in the elastic map and R2 distribution have been
marked with blank white spaces.
Fig. 5.17 shows the estimated Fung parameter b as a function of the estimated
indentation, computed using the inverse-FE routine. The inverse-FE routine was
invoked for correcting 268 force curves out of a total of 729 in the elastic map
shown in Fig. 5.16(a). Also shown in bold red is the predicted value bmax0.992 for a
linearly elastic response, as obtained from the sensitivity study. It can be seen
from Fig. 5.17(a) and Fig. 5.17(b) that hyperelastic force curves can be clearly
differentiated from the predicted linear elastic response (classification curve in bold
red). This would not have been possible without invocation of the inverse FE-
routine.
On closer examination of Fig. 5.17(b), it becomes apparent that there was a
significant population of force curves that showed very low values (∼ 0.001 − 0.1)
in the estimated parameter bˆ. One such representative force curve is shown in
Fig. 5.18. These force curves show surprisingly linear force-indentation curves, which
is probably due to material softening during the AFM loading. These force curves
did not fit well even after invocation of the inverse-FE routine.
147
0 200 400 600 800 1000
0
100
200
300
400
500
100 200 300 400 500 600 700 800
0
50
100
150
200
(a) (b)
Figure 5.17: (a) Variation of estimated Fung parameter bˆ from the Inverse FE
analysis as a function of the estimated indentation (using the contact estimate), (b)
zoomed image. Also overlaid in red is bmax0.992(∆).
0 100 200 300 400
0
20
40
60
80
100
120
Rsq:0.97703
Indentation (nm)
Fo
rc
e 
(n
N)
Data
Inverse FE
Dimitriadis
Figure 5.18: Representative force-indentation plot with low estimated b value (Eˆ0, bˆ)
= (153.21 kPa,0.005) with R2FE = 0.977. The Dimitriadis solution was E¯Dim was
146.14 kPa with R2Dim = 0.975.
148
5.8 Conclusions
AFM characterization studies on normal and cancerous breast tissue specimens
for histopathological inference necessitates very large scale AFM indentation data
acquisition and processing. In this chapter, some of the key issues in automating
AFM data analysis on spatially heterogeneous tissue specimens were addressed. The
major highlights of this chapter are enumerated as follows:
1. The applicability of the Hertz contact model in extracting elastic modulus from
AFM force curves on tissue specimens were highlighted and existing research
that overcame Hertzian limitations were reviewed.
2. To account for hyperelasticity beyond neo-Hookean approximation, the Fung
type phenomenological hyperelastic model was proposed to extract material
properties from highly nonlinear force curves.
3. Through a sensitivity study, bounds of the Fung model parameter b was ob-
tained that approximated a linear elastic tissue response.
4. A modified version of the BDP contact estimation algorithm[1] was proposed,
that could estimate the contact point accurately in nonlinear force curves.
5. An automation strategy was proposed that could accurately estimate material
properties and the contact point from simulated force curves with varying
degrees of hyperelasticity.
6. On implementation on an entire ROI, it was observed that the proposed au-
tomation strategy could significantly improve fit quality in the elastic map of
149
the ROI, thereby validating the research hypothesis of this chapter.
There are some limitations associated with the proposed approach discussed
in this chapter, namely:
1. In this analysis, a constant value of the indenter radius (R = 2.5µm) was
used. Varying the indenter size is likely to impact the results of the sensitivity
studies in section 5.5.3. However, the probe bead size is unlikely to affect the
major conclusions of this study, since it would still be possible to differentiate
between a linear and an hyperelastic tissue response based on the estimated
Fung parameter, bˆ from the inverse FE analysis.
2. It was observed that for force curves possessing very high degrees of hyperelas-
ticity, the contact estimation did not perform satisfactorily. This was primarily
because the force-indentation contact relationship used in the regression model
was not capable of explaining very high nonlinear tissue response. Establish-
ing a contact model for the Fung type hyperelastic materials remains a subject
of future work.
3. Anisotropy of the tissue response in the direction of indentation was not con-
sidered in this work. Currently, force curves that display anisotropy (shown in
Fig. 5.7) are isolated from the data analysis procedure after visually locating
them. An automated strategy using multi-changepoint estimation might be a
promising solution to identify and isolate these force curves in an unsupervised
manner.
150
Chapter 6
Conclusions and Future Work
6.1 Conclusions
This dissertation summarizes research directed towards performing accurate
mechanical characterization on breast tissue specimens with the intent of establish-
ing quantitative mechanical markers for cancer progression in the human breast.
Due to variability in the tissue elasticity across patients, it is difficult to conclu-
sively attribute mechanical signatures to the onset and progression of cancer in the
tissue specimens, however, this dissertation presents extensive work that has been
carried out towards achieving that goal. Several challenges were encountered while
carrying out preliminary AFM characterization experiments on tissue, which served
as the basis for chapter 2 of this dissertation. In the subsequent chapters 3-5, each
of these challenges were carefully examined and their solutions proposed. Based on
the research work done and experimental results obtained in this dissertation, the
following conclusions can be made:
1. Breast tissue specimens obtained from biopsy samples are highly heteroge-
neous in their composition and therefore, accurate mechanical characterization
should incorporate elastic map acquisitions (as discussed in chapter 3), instead
of random pointwise sampling in the annotated areas of the tissue specimen
(as discussed in chapter 2).
151
2. Obtaining representative tissue samples for quantifying mechanical property
alterations in cancerous tissue compared to normal tissue is a very challenging
exercise. As observed in the results in chapter 3, tissue elasticity is a strong
function of the individual patient from whom the biopsies were taken. It is
therefore recommended that for future characterization studies, normal and
cancer cores be extracted from the same patient.
3. Stromal tissue was consistently found to be more stiff compared to epithelial
regions in all the characterization experiments that were conducted. It is
therefore recommended that these regions need not be sampled as much as
epithelial regions for future characterization studies.
4. A probabilistic approach for tissue elastic modulus estimation provides a greater
degree of certainty in the characterization results since it incorporates errors
in the AFM probe calibration protocol and contact point uncertainty.
5. The direct implication of tissue heterogeneity is that AFM force curves ac-
quired during raster-scanning need to be carefully post-processed to identify
contribution of nonlinear elasticity in the force curves. Least square fitting of
the raw AFM data to existing contact models is attractive from the computa-
tional point of view; however, as observed in the results in chapter 5, a single
contact model does not account for all force curves obtained during raster
indentations on the tissue specimens.
152
6.2 Contributions
The contributions presented in this dissertation can be summarized as follows:
1. Development of a long range positioning system: A semi-automated position-
ing system has been proposed that overcomes some of the shortcomings of
using manual tissue-AFM probe alignment and manual registration methods.
Using image feedback, it was demonstrated that tissue ROIs beyond the range
of conventional AFM X-Y stages could be raster-scanned. It was also shown
that manual registration needed to be performed only twice during the AFM
experiments using the proposed system, in contrast to repetitive manual reg-
istration that needed to be periodically carried out to ensure that the probe
was localized correctly with the stained image.
2. EIV-Bayesian Changepoint algorithm: An integrated probabilistic approach
for estimating material properties from raw AFM data after incorporating
errors in contact estimation and AFM probe spring constant calibration vari-
ability was proposed. It was shown through extensive analysis that the algo-
rithm responded satisfactorily to varying degrees of calibration errors and the
accuracy of the results was validated on simulated AFM force curves.
3. Hyperelastic constitutive modeling of breast tissue: It is generally the practice
for AFM data analysis procedures to curve-fit bulk AFM data to analytical
contact models. While computationally inexpensive, might lead to inaccurate
estimates of the tissues’ elastic modulus when the tissue response is sharply
153
nonlinear. Using finite element modeling techniques, it was shown possible to
account for material nonlinearity through an exponential hyperelastic material
model, which significantly improved characterization accuracy.
6.3 Future Work
The possible directions for future work in this area are:
1. Characterization of stained tissue specimens: In this dissertation, AFM-based
tissue characterization was restricted to unstained tissue specimens to preclude
any tissue alterations due to staining. If stained tissue response parallels
those in unstained tissue, then stained tissue specimens could be used for
characterization. This would additionally aid in manual registration of the
annotated images with the tissue used for AFM probing.
2. Characterization of normal and cancerous specimens from the same patient:
Future work might be directed towards patient-specific characterization of
normal and cancerous tissue specimens.
3. Automated registration: Further throughput improvement with AFM charac-
terization could be achieved if the stained annotations could be automatically
registered with the tissue used for AFM probing. This can be performed using
multimodal registration techniques.
4. Non-parametric contact estimation: In the contact estimation strategies pro-
posed in chapters 4 and 5, an explicit contact model needed to be used. Con-
154
tact estimation strategies that do not require a contact model might be par-
ticularly attractive for estimating the contact point in force curves that are
not described by existing contact models.
5. Integrated computational framework for tissue characterization: In this dis-
sertation, the various steps involved in AFM characterization (probe spring
constant, contact estimation and tissue constitutive modeling) were treated in
a piecemeal manner. This was primarily because the EIV-Bayesian Change-
point algorithm (in chapter 4) required an explicit contact model and there
were no such models for highly nonlinear force curves, as seen in chapter 5. Fu-
ture work might be directed towards developing contact models for exponential
hyperelastic response on the lines of those proposed by Lin et al. [74] which
also accounts for finite-indentation. This would allow uncertainty quantifica-
tion of the tissue properties using the probabilistic tools developed in chapter
4. Confidence intervals thus generated could be used in interval-uncertainty
incorporated hypothesis tests, as suggested in [56] for inference purposes.
The experimental and computational strategies proposed in this dissertation
could have a significant impact on high-throughput quantitative studies of bioma-
terials, which could elucidate various disease mechanisms that are phenotyped by
their mechanical signatures.
155
Appendix A
Finite indentation response of Hertz and Dimitriadis contact models
In this appendix, the ability of Hertz and Dimitriadis contact models to ac-
count for finite indentation on thin specimens is investigated by comparing it with
a finite element (FE) response.
Simulated contact forces obtained by performing spherical indentation on a
linearly elastic FE model of a tissue specimen 4µm thick with varying elastic mod-
ulus E (from 10 kPa to 2500 kPa) were fitted with Dimitriadis and Hertz contact
models, resulting in the fitted moduli EDim and EHertz respectively. The normal-
ized parameters EDim/E and EHertz/E were obtained for varying indentation depth
(from ∆ = 50 nm to ∆ = 1400 nm). Two probe bead sizes R = 2.5µm and
R = 5.0µm were used for investigating the impact of indenter size on predictions
from Hertz and Dimitriadis contact models.
Fig. A.1 summarizes the results of the analysis performed. Figs. A.1(a) shows
the variation of EHertz/E and EDim/E for a probe bead radius of R = 2.5 µm
with increasing indentation depth and the corresponding fit quality (R2) is shown
in Fig. A.1(b). It becomes clear that the Hertz model significantly overestimates
the elastic modulus with increasing indentation depth due to the violation in in-
finitesimal assumptions. The corresponding fit quality too shows deteriorates with
increasing indentation. The overestimation in the Hertz response worsens further
156
for a larger bead size of R = 5.0 µm [Figs. A.1(c) and (d)].
0
500
1000
1500
0
1000
2000
1
1.5
2
2.5
3
 EDim/E
EHertz/E
0
500
1000
1500
0
1000
2000
1
1.5
2
2.5
3
0
500
1000
1500
0
1000
2000
3000
0.975
0.98
0.985
0.99
0.995
1
1.005
R2Dim
R2Hertz
0
500
1000
1500
0
1000
2000
3000
0.975
0.98
0.985
0.99
0.995
1
1.005
R2Dim
R2Hertz
EDim/E
EHertz/E
(c) (d)
(a) (b)
Figure A.1: (a) Variation of EHertz/E and EDim/E for R = 2.5 µm, (b) Corre-
sponding R2 fit quality, (c) Variation of EHertz/E and EDim/E for R = 5.0 µm and
(d) Corresponding R2 fit quality.
In contrast, it is clear that for both bead sizes, the Dimitriadis model accu-
rately estimates the elastic modulus (EDim/E ∼ 1) for shallow as well as deep in-
dentations. This analysis conclusively proves the applicability of Dimitriadis model
in estimating tissue elastic modulus from AFM indentation on thin specimens such
as the breast tissue specimens studied in this work.
157
Bibliography
[1] K.D. Costa, A.J. Smith, and F.C.P. Yin. Non-Hertzian Approach to Analyz-
ing Mechanical Properties of Endothelial Cells Probed by Atomic Force Mi-
croscopy. Journal of Biomechanical Engineering, 128(2):176–184, 2006.
[2] P.P. Rosen. Rosen’s Breast Pathology. Lippincott Williams & Wilkins, third
edition, 2008.
[3] Cancer facts & figures 2013. Atlanta: American Cancer Society, 2013.
[4] D.J. Foran, W. Chen, and L. Yang. Automated image interpretation and
computer-assisted diagnostics. Analytical Cellular Pathology, 34:279–300, 2011.
[5] J Ophir, S K Alam, B Garra, F Kallel, E Konofagou, T Krouskop, and T Vargh-
ese. Elastography: ultrasonic estimation and imaging of the elastic properties
of tissues. In In Proc. Instn. Mech. Engrs., Part H: Journal of Engineering in
Medicine, pages 203–233, 1999.
[6] A. Manduca, T.E. Oliphant, M.A. Dresner, J.L. Mahowald, S.A. Kruse, E. Am-
romin, J.P. Felmlee, J.F. Greenleaf, and R.L. Ehman. Magnetic resonance elas-
tography: non-invasive mapping of tissue elasticity. Medical Image Analysis,
5(4):237–254, 2001.
[7] A.J. Fredman, J.L. Frolik, and B.S. Garra. Lung strain profiles using com-
puted tomography elastography. In Annual Conference Proceedings of IEEE
Engineering in Medicine and Biology Society, pages 1545–1548, 2004.
[8] K.D. Costa. Single-cell elastography: Probing for disease with the atomic force
microscope. Disease Markers, 19(2):139–154, 2003.
[9] Q.S. Li, G.Y.H. Lee, C.N. Ong, and C.T. Lim. AFM indentation study of
breast cancer cells. Biochem. Biophys. Res. Commun., 374(4):609–613, 2008.
[10] M. Lekka, P. Laidler, D. Gil, J. Lekki, Z. Stachura, and A. Z. Hrynkiewicz.
Elasticity of normal and cancerous human bladder cells studied by scanning
force microscopy. European Biophyics Journal, 28(4):312–316, 1999.
[11] S. Park, D. Koch, R. Cardenas, J. Kas, and C.K. Shih. Cell motility and local
viscoelasticity of fibroblast. Biophysical Journal, 89:4330–4342, 2005.
[12] J. Katsantosis, A. Tosca, S.B. Koukouritaki, P.A. Theodoropoulos, A. Gravanis,
and C. Stournaras. Differences in the G/total actin ratio and microfilament
stability between normal and malignant keratinocytes. Cell Biochem Funct,
12:267–274, 1994.
158
[13] M.J. Paszek, N. Zahir, K.R. Johnson, A. Gefen, C.A. Reinhart-King, and S.S.
Boettiger. Tensional homeostasis and the malignant phenotype. Cancer Cell,
8(3):241–254, 2005.
[14] G.J. Mizejewski. Role of integrins in cancer: survey of expression patterns.
In Proceedings of the Society for Experimental Biology and Medicine, pages
124–138, 1999.
[15] D.T. Butcher, T. Alliston, and V.M. Weaver. A tense situation: forcing tumour
progression. Nature Rev. Cancer, 9(2):108–122, 2009.
[16] J. Kononen, L. Bubendorf, A. Kallionimeni, M. Brlund, P. Schraml,
S. Leighton, J. Torhorst, M.J. Mihatsch, G. Sauter, and O.P. Kallionimeni.
Tissue microarrays for high-throughput molecular profiling of tumor specimens.
Nature Medicine, 4:844–847, 1998.
[17] D. Kim, P.K. Wong, J. Park, A. Levchenko, and Y. Sun. Microengineered
platforms for cell mechanobiology. Annual Review of Biomedical Engineering,
11:202–233, 2009.
[18] A. Alessandrini and P. Facci. AFM : a versatile tool in biophysics. Measurement
Science and Technology, 16:R65–R92, 2005.
[19] A. Pillarisetti. Mechanical Characterization and Manipulation of Biological
Cells. PhD thesis, University of Maryland, College Park, 2008.
[20] W. Chen and D.J. Foran. Advances in tissue microarray technology: Towards
improved understanding and diagnosis of cancer. Anal. Chim. Acta, 564:74–81,
2006.
[21] W. Chen, D.J. Foran, and M. Reiss. Unsupervised imaging, registration and
archiving of tissue microarrays. In Proceedings of AMIA Symp, pages 136–139,
2002.
[22] W. Chen, D.J. Foran, and M. Reiss. A prototype for unsupervised analysis of
tissue microarrays for cancer research and diagnostics. IEEE Transactions on
Information Technology in Biomedicine, 8:89–96, 2004.
[23] L. Yang, O. Tuzel, W. Chen, P. Meer, G. Salaru, L.A. Goodell, and D.J. Foran.
Pathminer: A web-based tool for computer-assisted diagnostics in pathology.
IEEE Transactions on Information Technology in Biomedicine, 13(3):291–299,
2009.
[24] A. Pillarisetti, J.P. Desai, H. Ladjal, A. Schiffmacher, A. Ferreira, and C.L.
Keefer. Mechanical phenotyping of mouse embryonic stem cells: increase in
stiffness with differentiation. Cellular Reprogramming, 13(4):261–269, 2009.
[25] J.L. Hutter and J. Bechhoeffer. Calibration of atomic-force microscope tips.
Review Sci. Instrum., 64(7):1868–1873, 1993.
159
[26] K.L. Johnson. One hundred years of hertz contact. Proc Instn Mech Engrs,
196(1):363–378, 1982.
[27] J.E. Sader, I. Larson, P. Mulvaney, and L.R. White. Method for the calibration
of atomic force microscope cantilevers. Rev. Sci. Instrum., 66(7):3789–3798,
1995.
[28] R.S. Gates and J.R. Pratt. Accurate and precise calibration of AFM cantilever
spring constants using laser Doppler vibrometry. Nanotechnology, 23(37):1–12,
2012.
[29] David C. Lin, Emilios K. Dimitriadis, and Ferenc Horkay. Robust Strategies
for Automated AFM Force Curve Analysis-I. Non-adhesive Indentation of Soft,
Inhomogeneous Materials. Transactions of the ASME, 129(3):430–440, 2007.
[30] A.S. Naini, R.V. Patel, and A. Samani. Measurement of Lung Hyperelastic
Properties Using Iinverse Finite Element Approach. IEEE Transactions on
Biomedical Engineering, 58(10):2852–2859, 2011.
[31] H. Hagaa, S. Sasakia, K. Kawabataa, E. Itob, T. Ushikic, and T. Sambon-
gia. Elasticity mapping of living fibroblasts by AFM and immunofluorescence
observation of the cytoskeleton. Ultramicroscopy, 82(1-4):253–258, 2000.
[32] J.Shi and C. Tomasi. Good features to track. In Computer Vision and Pattern
Recognition, pages 593–600, 1994.
[33] D.G. Lowe. Distinctive Image Features from Scale-Invariant Keypoints. Inter-
national Journal of Computer Vision, 60(2):91–110, 2004.
[34] M.S. Nixon and A.S. Aguado. Feature Extraction & Image Processing. Aca-
demic Press, second edition, 2008.
[35] Gary Bradski and Adrian Kaehler. Learning OpenCV. 2008.
[36] Y. Sun and B. Nelson. Biological Cell Injection Using an Autonomous Mi-
croRobotic System. The International Journal of Robotics Research, 21(10-
11):861–868, 2002.
[37] F. Liu and D.J. Tschumperlin. Micro-mechanical characterization of lung tissue
using atomic force microscopy. Journal of Visualized Experiments, 54(2911),
2011.
[38] B. Capella, P. Baschieri, C. Frediani, P. Miccoli, and C. Ascoli. Force-Distance
Curves by AFM. IEEE Engineering in Medicine and Biology, 16(2):58–65,
1997.
[39] M.F. Jaasma, W.M. Jackson, and T.M. Keaveny. Measurement and Character-
ization of Whole-Cell Mechanical Behavior. Annals of Biomedical Engineering,
34(5):748–758, 2006.
160
[40] L.R. Nyland and D.W. Maughan. Morphology and transverse stiffness of
drosophila myofibrils measured by atomic force microscopy. Biophysical Jour-
nal, 78(3):14901497, 2000.
[41] M. Radermacher. Measuring the elastic properties of living cells by the atomic
force microscope. Methods Cell Biol, 68:67–90, 2002.
[42] G.T. Charras, P.P. Lehenkari, and M.A. Horton. Atomic force microscopy
can be used to mechanically stimulate osteoblasts and evaluate cellular strain
distributions. Ultramicroscopy, 86:85–95, 2001.
[43] M. Radmacher, M. Fritz, C. Kacher, J. Cleveland, and P. Hansma. Measuring
the Viscoelastic Properties of Human Platelets with the Atomic Force Micro-
scope. Biophysical Journal, 70(1):556–567, 1996.
[44] E.K. Dimitriadis, F. Horkay, J. Maresca, B. Kachar, and R.S. Chadwick. De-
termination of elastic moduli of thin layers of soft material using the atomic
force microscope. Biophysical Journal, 82(5):2798–2810, 2002.
[45] S.L. Crick and F.C.-P. Yin. Assessing Micromechanical Properties of Cells with
Atomic Force Microscopy: Importance of the Contact Point. Biomechanics and
modeling in Mechanobiology, 6(3):199–210, 2007.
[46] Pavel Polyakov, Charles Soussen, Junbo Duan, Jerome F. L. Duval, David Brie,
and Gregory Francius. Automated force volume image processing for biological
samples. PLoS ONE, 6(4):e18887, 2001.
[47] E.A-Hassan, W.F. Heinz, M.D. Antonik, N.P. D’Costa, S. Nageswaran, C.A.
Schoenenberger, and J.H. Hoh. Relative microelastic mapping of living cells by
atomic force microscopy. Biophysical Journal, 74(3):1564–1578, 1998.
[48] Daniel Rudoy, Shelten G. Yuen, Robert D. Howe, and Patrick J. Wolfe.
Bayesian change-point analysis for atomic force microscopy and soft material
indentation. Journal of the Royal Statistical Society: Series C, 59(4):573–593,
2010.
[49] J.P. Cleveland, S. Manne, D. Bocek, and P.K. Hansma. A nondestructive
method for determining the spring constant of cantilevers for scanning force
microscopy. Review of Scientific Instruments, 64:403, 1993.
[50] A. Torii, M. Sasaki, K. Hane, and S. Okuma. A method for determining the
spring constant of cantilevers for atomic force microscopy. Measurement Science
and Technology, 7(2):179–184, 1996.
[51] R. Proksch, T. Schaffer, J.P. Cleveland, R.C Callahan, and M.B. Viani. Finite
optical spot size and position corrections in thermal spring constant calibration.
Nanotechnology, 15(9):1344–1350, 2004.
161
[52] N.A. Burnham, X. Chen, C.S. Hodges, G.A. Mateil, E.J. Thoreson, C.J.
Roberts, M.C. Davies, and S.J.B. Tendler. Comparison of calibration meth-
ods for atomic-force microscopy cantilevers. Nanotechnology, 14(1):1–6, 2003.
[53] T. Pirzer and T. Hugel. Atomic force microscopy spring constant determination
in viscous liquids. Review of Scientific Instruments, 80(035110):1–6, 2009.
[54] B. Ohler. Cantilever constant calibration using laser Doppler vibrometry. Re-
view of Scientific Instruments, 78(6):063701, 2007.
[55] D. Kiracofe and A. Raman. Quantitative force and dissipation measurements
in liquids using piezo-excited atomic force microscopy: a unifying theory. Nan-
otechnology, 22:1–11, 2011.
[56] V. Kreinovich, H.T. Nguyen, and S. Niwitpong. Statistical hypothesis testing
under interval uncertainty: An overview. International Journal of Intelligent
Technology and Applied Statistics, 1(1):1–32, 2008.
[57] R. Roy, W. Chen, L. Cong, D.J. Foran, and J.P. Desai. Probabilistic Estimation
of Mechanical Properties of Biomaterials Using Atomic Force Microscopy. IEEE
Tranactions on Biomedical Engineering, 2013.
[58] R. Wagner, R. Moon J. Pratt G. Shaw, and A. Raman. Uncertainty quan-
tification in nanomechanical measurements using the atomic force microscope.
Nanotechnology, 22(45):1–9, 2011.
[59] R. Long, M. Hall, M. Wu, and C-Y. Hui. Effects of Gel Thickness on Microscopic
Indentation Measurements of Gel Modulus. Biophysical Journal, 101(3):642–
650, 2011.
[60] C-L. Cheng and J.W. Van Ness. Statistical Regression with Meaurement Error
: Kendall’s Library of Statistics 6. Wiley, first edition, 1999.
[61] B.P. Carlin, A.E. Gelfand, and A.F.M. Smith. Hierarchical Bayesian Analysis
of Changepoint Problems. Journal of the Royal Statistical Society. Series C
(Applied Statistics), 41(2):389–405, 1992.
[62] A. Gelman, J.B. Carlin, H.S. Stern, and D.B. Rubin. Bayesian Data Analysis.
Chapman and Hall, first edition, 1995.
[63] D. Stephens. Bayesian Retrospective Multiple-Changepoint Identification.
Journal of the Royal Statistical Society, 43(1):159–178, 1994.
[64] C. Gossl and H.Kuchenoff. Bayesian analysis of logistic regression with an
unknown change point and covariate measurment error. Statistics in Medicine,
20(3):3109–3121, 2001.
162
[65] R. Roy, W. Chen, L. Cong, L.A. Goodell, D.J. Foran, and J.P. Desai.
A Semi-Automated Positioning System for contact-mode Atomic Force Mi-
croscopy (AFM). IEEE Transactions on Automation Science and Engineering,
10(2):462–465, 2013.
[66] K.D. Costa and F.C.P. Yin. Analysis of indentation: Implications for measuring
mechanical properties with atomic force microscopy. Journal of Biomechanical
Engineering, 121(5):462–471, 1999.
[67] L. Sirghi, J. Ponti, F. Broggi, and F. Rossi. Probing elasticity and adhesion of
live cells by atomic force microscopy indentation. European Biophysics Journal,
37(6):935–945, 2008.
[68] David C. Lin, Emilios K. Dimitriadis, and Ferenc Horkay. Robust Strategies
for Automated AFM Force Curve Analysis-II - Adhesion-Inuenced Indentation
of Soft, Elastic Materials. Transactions of the ASME, 129(6):904–912, 2007.
[69] R.E. Mahaffy, C.K. Shih, F.C. MacKintosh, and J. Kas. Scanning Probe-
Based Frequency-Dependent Microrheology of Polymer Gels and Biological
Cells. Physical Review Letters, 85(4):880–883, 2000.
[70] J. Domke and M. Radmacher. Measuring the Elastic Properties of Thin Poly-
mer Films with the Atomic Force Microscope. Langmuir, 14:3320–3325, 1998.
[71] K. Sangpradit, H. Liu, L. Seneviratne, and K. Althoefer. Tissue Identifica-
tion using Inverse Finite Element Analysis of Rolling Indentation. In In IEEE
International Conference on Robotics and Automation, pages 1250–1255, 2009.
[72] S. Park, K.D. Costa, and G.A. Ateshian. Microscale Frictional Response of
Bovine Articular Cartilage from Atomic Force Microscopy. Journal of Biome-
chanics, 37(11):1679–1687, 2004.
[73] J.A.C. Santos, L.M.Rebelo, A.C. Araujo, E.B. Barrosa, and J.S. de Sousa.
Thickness-corrected model for nanoindentation of thin films with conical in-
denters. Soft Matter, 8(16):4441–4448, 2012.
[74] D.C. Lin, D.I. Shreiber, E.K. Dimitriadis, and F. Horkay. Spherical indentation
of soft matter beyond the Hertzian regime: numerical and experimental vali-
dation of hyperelastic models. Biomech Model Mechanobiol, 8:345–358, 2009.
[75] D.C. Lin, D.I. Shreiber, E.K. Dimitriadis, and F. Horkay. Elasticity of rubber-
like materials measured by AFM nanoindentation. eXPRESS Polymer Letters,
1(9):576–584, 2007.
[76] D. Tabor. The Hardness of Materials. Clarendon Press, Oxford, first edition,
1951.
163
[77] B. Oommen and K.J. Van Vilet. Effect of nanoscale thickness and elastic
nonlinearity on measured properties of polymeric films. Thin Solid Films, 513(1-
2):235–242, 2006.
[78] H. Ladjal, J-Luc Hanus, A. Pillarisetti, C. Keefer, A. Ferreira, and J.P. De-
sai. Atomic Force Microscopy-Based Single-Cell Indentation: Experimentation
and Finite Element Simulation. In IEEE/RSJ International Conference on
Intelligent Robots and Systems, pages 1326–1332, 2009.
[79] M.G. Zhong, Y.P. Cao, G.Y. Li, and X.Q. Feng. Spherical indentation method
for determining constitutive paramters of hyperelastic soft materials. Biomech
Model Mechanobiol, 2013.
[80] M.C. Frassanito, L. Lamberti, A. Boccaccio, and C. Pappalettere. Discussion
on hybrid approach to determination of cell elastic properties. In Optical Mea-
surements, Modeling, and Metrology, Volume 5, Conference Proceedings of the
Society for Experimental Mechanics Series, pages 119–124, 2011.
[81] P.J. Thurner. Atomic force microscopy and indentation force measurement of
bone. Nanomedicine and Nanobiotechnology, 1(6):624–649, 2009.
[82] T. Yamada, Y. Kunioka, J. Wakayama, M. Aimi, Y.S. Noguchi amd
N. Akiyama, and T. Kayamori. Molecular organizations of myofibrils of skeletal
muscle studied by atomic force microscopy. Advances in experimental medicine
and biology, 538:285–294, 2003.
[83] H.K. Herisa, A.K. Miria, U. Tripathy, F. Barthelata, and L. Mongeaua. Inden-
tation of poroviscoelastic vocal fold tissue using an atomic force microscope.
Journal of the Mechanical Behavior of Biomedical Materials, 28:383–392, 2013.
[84] G.U. Unnikrishnan, V.U. Unnikrishnan, and J.N. Reddy. Constituitive material
modeling of cell: A micromechanics approach. ASME Journal of Biomechanical
Engineering, 129(3):315–323, 2007.
[85] J.A. Nelder and R. Mead. A simplex method for function minimization. Com-
puter Journal, 7(4):308–313, 1965.
[86] O.H. Yeoh. Some forms of the strain energy function for rubber. Rubber Chem-
istry and Technology, 66(5):754–771, 1993.
[87] Y.C. Fung. Biomechanics: Mechanical Properties of Living Tissues. Springer,
second edition, 1993.
164