Fitness AFM Example
From QuB
In this article we use Fitness curve fitting software to analyze force-distance curves from an atomic force microscope (AFM) probing a living cell. We will derive a custom fitting function and use it to simultaneously fit the contact point (cell height) and the cell stiffness (elasticity).
Contents |
Data Preparation
Each trace was saved in a tab-separated text file with descriptive column headers:
Time Force Height 103.218260 -43.41 -2.695 103.218860 18.96 -2.69 103.219460 1.483 -2.685 103.220060 15.19 -2.68 103.220660 -3.315 -2.675
All of the trace files from one experiment were kept together in one folder.
Here is a trace loaded into Fitness (File -> Open data file...):
Notice how the data from 50% to 100% is selected for fitting, and "Force" and "Height" are chosen for the X and Y axes.
Custom Fit Curve
As described in (Spagnoli et al, Physical Review E 2008), the Hertz model of a hard sphere indenting an elastic solid surface "predicts that
,
where R is the radius of the indenter, z is the indentation distance relative to the plane surface, E is Young’s modulus, and ν is Poisson’s constant. We assumed </math>\nu</math> = 0.3 and using a nonlinear optimizer (QFC aka Fitness) fit each trace to the Hertz response with a linear baseline to allow for possible drift in a given F-z curve (parameters m and b):
F = mz + b + k * (z − z0)1.5 * H(z − z0),
where H is the Heaviside step function and 
."
To enter this curve in Fitness, we translate it to a Python expression in terms of x (not z):
m*x + b + H(x-x0)*k*(x-x0)**1.5
and replace the Heaviside function with its Python "short-circuit" equivalent:
m*x + b + (x>=x0 and k*(x-x0)**1.5 or 0)
We enter this directly into the "Curve" window in Fitness. Then we initialize the parameters m, b, x0 and k to reasonable values, check the box next to each parameter name, and click "Fit."
Finally, we press "Keep It" to add the parameters and errors to the "Fit Table"
To save the formula for later use, click "Presets..." and "Add to menu..."
Discussion
In the data files, "x" is mapped to probe advance towards the substrate surface and "y" is mapped to force. x0 is the contact point. The first part (mx+b) simply calculates a straight line fit to a sloping baseline. The second part needs a bit of clarification, especially if you are not familiar with python. The use of the arithmetic "and/or" in python is specifically designed to return the first valid non-zero term as the arithmetic result of the operator. The rest of the expression will be ignored after that point. Thus, when the probe is higher than the contact point, x is less than x0, the "and" expression is false and the "or" part of the expression (zero) is returned. As soon as x reaches x0, x>=x0 becomes true and the value of k*(x-x0)**1.5 is returned and added to the baseline slope. Note that the odd case where x exactly equals x0 works out just fine since the "or" clause will return 0 in this case, which is the result of the k*(x-x0)**1.5 calculation anyway.
The least squares approximation then estimates m, b, x0 and k. Young's modulus was then calculated from k using E=3*(1-nu**2)/(4*sqrt(R)), with nu=0.36 and R=7.5 microns, the diameter of the glass bead which was glued to the tip of the cantilever.
At first examination, it is hard to believe that this mess works at all, given the function switch in the middle of the residuals evaluation. However, even though the second derivative is discontinuous, the first derivative is not, so estimation of the next value of the parameter succeeds quite nicely. The major limitation is that the value of x0 must be constrained to lie on the interval of the range of x values (actually it only needs to be lower that the largest x value). As soon as the value of x0 gets to be greater than the largest x value, then evaluation of the k*(x-x0)**1.5 ceases and the the derivative of the function with respect to x0 becomes 0 so the fit can never bring x0 back into range. The result will simply be the best fit straight line. However, as long as x0<xmax, then the fit generally succeeds. This was only an issue for us when large numbers of curves were being fit and/or when xmax shortens significantly from curve to curve. Our solution to this problem was to always initialize x0 to a value at about 90% of the x-range. Sometimes it was necessary to also use the Strategy feature to fit twice: first with x0 constrained to this fixed value, then enable x0 and fit the second time.
I've removed the cosine term since I suspect that the baseline vibration problem that we experienced was atypical. Especially since we often had days when the data was very well represented using a simple straight line slope in the baseline leading up to contact. Also, to simplify the example, I've removed the part of the description that recalculates the error weighting terms, since these were only needed when I needed to fit several hundred curves in an automatic fashion. I can pass on the details of how we did the weights if you want them. In most cases, simply defining a fixed weighting function works just fine and even that is not needed most of the time.
Multi-file Fitting
The "Fit All" button fits the curve to each file in the current folder. To start, we open one of the files and fit it, to initialize the curve and parameters. Then press "Fit All." When it's done we use the "Fit Table" window to browse the results.
To fit all traces correctly, we had to change the Fit button's "Strategy" to fit twice: first with x0 constrained to this fixed value, then enable x0 and fit the second time.
- under Curve, click "Strategy"
- click "Add a step" to repeat the default Fit process twice
- in the first step, remove 'x0' from the list of 'variables'
- in the first step, replace "x0=Initial" with "x0=0.0"
- check the box next to "Strategy" to turn it on
- click "Fit All"
The strategy shown below initializes x0 to a value at 90% of the x range:
