docs: document surface calibration

changelog.md

@@ -19,6 +19,7 @@
 * Plot the active calibration peak for a calibration result using `show_active_peak=True` with [`CalibrationResults.plot()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults.html#lumicks.pylake.force_calibration.power_spectrum_calibration.CalibrationResults.plot).
 * Added function to import a [`KymoTrackGroup`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.kymotracker.kymotrack.KymoTrackGroup.html) from a `CSV` file using [`load_tracks`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.load_tracks.html).
 * Added function to load tracks into the kymotracker widget using [`KymoWidgetGreedy.load_tracks()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.KymoWidgetGreedy.html#lumicks.pylake.KymoWidgetGreedy.load_tracks).
+* Added [`lk.touchdown()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.touchdown.html#lumicks.pylake.touchdown) to find the height above the surface using the axial force signal and added an example notebook demonstrating surface calibration.
 
 #### Improvements
 

docs/api.rst

@@ -42,6 +42,7 @@ Force calibration
     force_calibration.power_spectrum.PowerSpectrum
     force_calibration.power_spectrum_calibration.CalibrationResults
     force_calibration.calibration_models.DiodeCalibrationModel
+    force_calibration.touchdown.TouchdownResult
 
     :template: function.rst
 
@@ -50,6 +51,7 @@ Force calibration
     fit_power_spectrum
     viscosity_of_water
     density_of_water
+    touchdown
     coupling_correction_2d
 
 

docs/examples/index.rst

@@ -22,5 +22,6 @@ For all of the examples, it is assumed that the following lines precede any othe
     cas9_kymotracking/cas9_kymotracking
     hairpin_fitting/hairpin_unfolding
     droplet_fusion/droplet_fusion
+    surface_calibration/surface_calibration.rst
     bead_coupling/coupling
     binding_lifetime/binding_lifetime

docs/refs.bib

@@ -32,6 +32,28 @@ @article{Kaczmarczyk2022
   doi = {10.1038/s41467-022-33503-6}
 }
 
+@article{neuman2004optical,
+  title={Optical trapping},
+  author={Neuman, Keir C and Block, Steven M},
+  journal={Review of scientific instruments},
+  volume={75},
+  number={9},
+  pages={2787--2809},
+  year={2004},
+  publisher={American Institute of Physics}
+}
+
+@article{neuman2005measurement,
+  title={Measurement of the effective focal shift in an optical trap},
+  author={Neuman, Keir C and Abbondanzieri, Elio A and Block, Steven M},
+  journal={Optics letters},
+  volume={30},
+  number={11},
+  pages={1318--1320},
+  year={2005},
+  publisher={Optica Publishing Group}
+}
+
 @article{Kochaniak2009,
   author = {Kochaniak, A B and Habuchi, S and Loparo, J J and Chang D J and Cimprich K A and Walter J C and van Oijen A M},
   title = {Proliferating Cell Nuclear Antigen Uses Two Distinct Modes to Move along DNA},

lumicks/pylake/__init__.py

@@ -25,6 +25,7 @@
 from .kymotracker.kymotracker import *
 from .piezo_tracking.baseline import *
 from .nb_widgets.range_selector import FdRangeSelector, FdDistanceRangeSelector
+from .force_calibration.touchdown import touchdown
 from .force_calibration.convenience import calibrate_force
 from .piezo_tracking.piezo_tracking import *
 from .nb_widgets.kymotracker_widgets import KymoWidgetGreedy
@@ -34,6 +35,7 @@
     density_of_water,
     viscosity_of_water,
     coupling_correction_2d,
+    surface_drag_correction,
 )
 from .force_calibration.power_spectrum_calibration import (
     fit_power_spectrum,

lumicks/pylake/force_calibration/calibration_models.py

@@ -222,6 +222,38 @@ def viscosity_of_water(temperature, molarity_nacl=None, pressure=None):
         return _poly((temperature + 273.15) / 300, bi, ai) * 1e-6
 
 
+def surface_drag_correction(distance_to_surface, bead_diameter, axial=False):
+    """Calculate a correction factor for the drag coefficient due to a nearby surface [1]_.
+
+    Parameters
+    ----------
+    distance_to_surface : array_like | float
+        Distance from the center of the bead to the surface
+    bead_diameter : float
+        Bead diameter
+    axial : bool
+        Compute the correction factor for axial drag
+
+    References
+    ----------
+    .. [1] Schäffer, E., Nørrelykke, S. F., & Howard, J. "Surface forces and drag coefficients of
+           microspheres near a plane surface measured with optical tweezers." Langmuir, 23(7), 3654-3665
+           (2007).
+    """
+    distance_to_surface = np.asarray(distance_to_surface)
+
+    if np.any(distance_to_surface <= 0.5 * bead_diameter):
+        raise ValueError(
+            f"Bead is inside the surface. you specified a distance of {distance_to_surface}, while"
+            f"the bead radius is {bead_diameter / 2}"
+        )
+
+    if axial:
+        return brenner_axial(distance_to_surface, bead_diameter / 2)
+    else:
+        return faxen_factor(distance_to_surface, bead_diameter / 2)
+
+
 def coupling_correction_2d(dx, dy, bead_diameter, is_y_oscillation=False, allow_rotation=True):
     """Calculates the coupling correction factor for a 2D problem.
 

lumicks/pylake/force_calibration/detail/drag_models.py

@@ -14,7 +14,7 @@ def faxen_factor(distance_to_surface_m, radius_m):
 
     Parameters
     ----------
-    distance_to_surface_m : float
+    distance_to_surface_m : array_like | float
         Distance from the center of the bead to the surface [m]
     radius_m : float
         Radius of the bead [m]
@@ -31,7 +31,7 @@ def faxen_factor(distance_to_surface_m, radius_m):
 
 
 def brenner_axial(distance_to_surface_m, radius_m):
-    """Brenner factor for lateral drag coefficient.
+    """Brenner factor for axial drag coefficient.
 
     This factor provides a correction to the drag force for a nearby wall.
 
@@ -41,7 +41,7 @@ def brenner_axial(distance_to_surface_m, radius_m):
 
     Parameters
     ----------
-    distance_to_surface_m : float
+    distance_to_surface_m : array_like | float
         Distance from the center of the bead to the surface [m]
     radius_m : float
         Radius of the bead [m]

lumicks/pylake/force_calibration/tests/test_hydro.py

@@ -3,6 +3,7 @@
 from lumicks.pylake.force_calibration.calibration_models import (
     ActiveCalibrationModel,
     PassiveCalibrationModel,
+    surface_drag_correction,
 )
 from lumicks.pylake.force_calibration.detail.drag_models import *
 from lumicks.pylake.force_calibration.detail.power_models import g_diode
@@ -492,6 +493,11 @@ def test_faxen(bead_radius, result):
     np.testing.assert_allclose(
         faxen_factor(bead_radius + np.arange(40, 250, 80) * 1e-9, bead_radius), result
     )
+    np.testing.assert_allclose(
+        surface_drag_correction(bead_radius + np.arange(40, 250, 80) * 1e-9, bead_radius * 2),
+        result,
+        False,
+    )
 
 
 @pytest.mark.parametrize(
@@ -505,6 +511,11 @@ def test_brenner(bead_radius, result):
     np.testing.assert_allclose(
         brenner_axial(bead_radius + np.arange(40, 250, 80) * 1e-9, bead_radius), result
     )
+    np.testing.assert_allclose(
+        surface_drag_correction(bead_radius + np.arange(40, 250, 80) * 1e-9, bead_radius * 2),
+        result,
+        True,
+    )
 
 
 def test_near_surface_consistency():

lumicks/pylake/force_calibration/touchdown.py

@@ -208,15 +208,23 @@ def cost(freq):
 
 @dataclass
 class TouchdownResult:
-    surface_position: float
-    focal_shift: float
+    surface_position: float | None
+    """Position where the coverslip surface touches the bead."""
+    focal_shift: float | None
+    """Focal shift of the optical setup."""
     nanostage_position: np.ndarray
+    """Nanostage position in microns."""
     axial_force: np.ndarray
+    """Axial force in pN."""
     surface_fit: np.ndarray
+    """Fit to the surface approach curve."""
     interference_nanostage: np.ndarray
+    """Fit to the interference pattern."""
     interference_force: np.ndarray
+    """Interference pattern."""
 
     def plot(self, legend=True):
+        """Plot the results of the touchdown analysis."""
         import matplotlib.pyplot as plt
 
         plt.plot(self.nanostage_position, self.axial_force, label="Axial force")
@@ -231,11 +239,54 @@ def plot(self, legend=True):
             self.surface_fit,
             label=f"Piecewise linear fit{'' if self.surface_position else ' (failed)'}",
         )
+        plt.xlabel("Nanostage position (µm)")
+        plt.ylabel("Axial force (pN)")
         if self.surface_position:
-            plt.axvline(self.surface_position, label="Determined surface position")
+            plt.axvline(
+                self.surface_position,
+                label="Determined contact position",
+                color="k",
+                linestyle="--",
+            )
         if legend:
             plt.legend()
 
+    def calculate_height(self, nanostage_position, bead_diameter):
+        r"""Calculates the height above the surface for a nanostage position and bead diameter.
+
+        Note that this returns the distance from the bead center to the surface, which is the
+        distance that you should supply to other force calibration routines.
+
+        Computes the following equation:
+
+        .. math::
+
+            d / 2 - \alpha_{shift} \left(z_{nanostage} - z_{surface}\right)
+
+        Here :math:`d` is the bead diameter, :math:`\alpha_{shift}` is the focal shift factor,
+        :math:`z_{nanostage}` is the nanostage position and :math:`z_{surface}` is the nanostage
+        position at which the bead and flowcell touch (surface-to-surface).
+
+        .. note ::
+
+            This function assumes that the trap z-position has not been changed after the
+            touchdown analysis and does not take into account errors in the estimate of the bead
+            radius.
+
+        Parameters
+        ----------
+        nanostage_position : array_like
+            Nanostage position at which we want to know the height above the surface.
+        bead_diameter : float
+            Bead diameter
+        """
+        if not self.surface_position or not self.focal_shift:
+            raise RuntimeError("Cannot compute height above the surface because the fit failed")
+
+        return 0.5 * bead_diameter + self.focal_shift * (
+            self.surface_position - np.asarray(nanostage_position)
+        )
+
 
 def touchdown(
     nanostage,
@@ -301,6 +352,7 @@ def touchdown(
         raise ValueError("Insufficient data available to fit touchdown curve")
 
     stage_trimmed, force_trimmed = nanostage[mask], axial_force[mask]
+    # Wavelength of the modulation pattern
     expected_wavelength = wavelength_nm * 1e-3 / 2 / refractive_index_medium
     bounds = np.array([0.7, 1.0001]) / expected_wavelength
     par, simulation = fit_damped_sine_with_polynomial(