Skip to content

Commit

Permalink
hm
Browse files Browse the repository at this point in the history
  • Loading branch information
JoepVanlier committed Oct 10, 2024
1 parent 6185d34 commit 7a26ab6
Show file tree
Hide file tree
Showing 9 changed files with 218 additions and 208 deletions.
2 changes: 2 additions & 0 deletions docs/theory/force_calibration/active.rst
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
.. _active_calibration_theory:

Active Calibration
------------------

Expand Down
4 changes: 4 additions & 0 deletions docs/theory/force_calibration/hyco.rst
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,10 @@ We can omit this effect by passing `fast_sensor=True` to the calibration models
Note however, that this makes using the hydrodynamically correct model critical, as the simple model doesn't actually capture the data very well.
The following example data acquired on a fast sensor will illustrate why::

filenames = lk.download_from_doi("10.5281/zenodo.7729823", "test_data")
f = lk.File("test_data/fast_measurement_25.h5")

# Decalibrate the force data
volts = f.force2y / f.force2y.calibration[0].force_sensitivity

shared_parameters = {
Expand Down
184 changes: 0 additions & 184 deletions docs/tutorial/force_calibration.rst
Original file line number Diff line number Diff line change
@@ -1,136 +1,3 @@
Calibration near the surface
----------------------------

So far, we have only considered experiments performed deep in bulk.
In reality, proximity of the flowcell surface to the bead leads to an increase in the drag force on the bead.
Such surface effects can be taken into account by specifying a distance to the surface (in microns)::

force_model = lk.PassiveCalibrationModel(
bead_diameter,
hydrodynamically_correct=True,
rho_sample=999,
rho_bead=1060.0,
distance_to_surface=5
)

As we approach the surface, the drag effect becomes stronger.
The hydrodynamically correct model is only valid up to a certain distance to the surface.
Moving closer, the frequency dependent effects become smaller than the overall drag effect and better models to approximate the local drag exist.
We can use this model by setting `hydrodynamically_correct` to `False`, while still providing a distance to the surface::

force_model = lk.PassiveCalibrationModel(
bead_diameter, hydrodynamically_correct=False, distance_to_surface=5
)

To summarize, the workflow can be visualized as follows:

.. image:: figures/force_calibration/surface_calibration_workflow.png

.. note::

Note that the drag coefficient `gamma_ex` that Pylake returns always corresponds to the drag coefficient extrapolated back to its bulk value.
This ensures that drag coefficients can be compared and carried over between experiments performed at different heights.
The field `local_drag_coefficient` contains an estimate of the local drag coefficient (at the provided height).

Axial Force
-----------

No hydrodynamically correct model is available for axial calibration.
However, models do exist for the dependence of the drag force on the distance to the surface.
Axial force calibration can be performed by specifying `axial=True`::

force_model = lk.PassiveCalibrationModel(bead_diameter, distance_to_surface=5, axial=True)

Active calibration with a single bead
-------------------------------------

Active calibration has a few benefits.
When performing passive calibration, we base our calculations on a theoretical drag coefficient which depends on parameters that are only known with limited precision:

- The diameter of the bead :math:`d` in microns.
- The dynamic viscosity :math:`\eta` in Pascal seconds.
- The distance to the surface :math:`h` in microns.

The viscosity in turn depends strongly on the local temperature around the bead, which is typically poorly known.

In active calibration, we oscillate the stage with a known frequency and amplitude.
This introduces an extra peak in the power spectrum which allows the trap to be calibrated without the assumptions of the theoretical drag coefficient.

Using Pylake, the procedure to use active calibration is not very different from passive calibration.
However, it does require some additional data channels as inputs.
In the next section, the aim is to calibrate the x-axis of trap 1.
We will consider that the nanostage was used as driving input.
Let's analyze some active calibration data acquired near a surface.
To do this, load a new file::

f = lk.File("test_data/near_surface_active_calibration.h5")
volts = decalibrate(f.force1x)
bead_diameter = f.force1x.calibration[0]["Bead diameter (um)"]
# Calibration performed at 1.04 * bead_diameter
distance_to_surface = 1.04 * bead_diameter

First we need to extract the nanostage data which is used to determine the driving amplitude and frequency::

driving_data = f["Nanostage position"]["X"]

For data acquired with active calibration in Bluelake, this will be a sinusoidal oscillation.
If there are unexplained issues with the calibration, it is a good idea to plot the driving signal and verify that the motion looks like a clean sinusoid::

plt.figure()
driving_data.plot()
plt.xlim(0, 0.1)
plt.ylabel(r"Nanostage position ($\mu$m)")
plt.show()

.. image:: figures/force_calibration/nanostage_position.png

Instead of using the :class:`~lumicks.pylake.PassiveCalibrationModel` presented in the previous section, we now use the :class:`~lumicks.pylake.ActiveCalibrationModel`.
We also need to provide the sample rate at which the data was acquired, and a rough guess for the driving frequency.
Pylake will find an accurate estimate of the driving frequency based on this initial estimate (provided that it is close enough)::

active_model = lk.ActiveCalibrationModel(
driving_data.data,
volts.data,
driving_data.sample_rate,
bead_diameter,
driving_frequency_guess=38,
distance_to_surface=distance_to_surface
)

We can check the determined driving frequency with::

>>> active_model.driving_frequency
38.15193077664577

Let's have a look to see if this peak indeed appears in our power spectrum.
To see it clearly, we reduce the blocking amount and show the spectrum all the way up to a frequency of 10 Hz::

show_peak = lk.calculate_power_spectrum(
volts.data, sample_rate=volts.sample_rate, num_points_per_block=5, fit_range=(10, 23000)
)

plt.figure()
show_peak.plot()
plt.show()

.. image:: figures/force_calibration/calibration_peak.png

The driving peak is clearly visible in the spectrum.
Next let's calculate the power spectrum we'll use for fitting.
It is important to *not* include the driving peak when doing this (the default will only go up to 100 Hz)::

power_spectrum = lk.calculate_power_spectrum(volts.data, sample_rate=volts.sample_rate)

We can now use this to fit our data::

calibration = lk.fit_power_spectrum(power_spectrum, active_model)
calibration

.. image:: figures/force_calibration/calibration_item_active.png

Analogous to the passive calibration procedure, we can specify `hydrodynamically_correct=True` if we wish to use the hydrodynamically correct theory here.
Especially for bigger beads this is highly recommended (more on this later).

Comparing different types of calibration
----------------------------------------

Expand Down Expand Up @@ -178,57 +45,6 @@ However, if we do not provide the height above the surface, we can see that the
Consequently, when this model is selected, this parameter affects both passive and active calibration.
For more information on this see the :doc:`theory section on force calibration</theory/force_calibration/force_calibration>` section.

.. _bead_bead_tutorial:

Active calibration with two beads far away from the surface
-----------------------------------------------------------

.. warning::

The implementation of the coupling correction models is still alpha functionality.
While usable, this has not yet been tested in a large number of different scenarios.
The API can still be subject to change *without any prior deprecation notice*!
If you use this functionality keep a close eye on the changelog for any changes that may affect your analysis.

When performing active calibration with two beads, we get a lower fluid velocity around the beads than we would with a single bead.
This leads to a smaller voltage readout than expected and therefore a higher displacement sensitivity (microns per volt).
Failing to take this into account results in a bias.
Pylake offers a function to calculate a correction factor to account for the lower velocity around the bead.

.. note::

For more information on how these factors are derived, please refer to the :ref:`theory<bead_bead_theory>` section on this topic.

Appropriate correction factors for oscillation in x can be calculated as follows::

factor = lk.coupling_correction_2d(dx=5.0, dy=0, bead_diameter=bead_diameter, is_y_oscillation=False)

Here `dx` and `dy` represent the horizontal and vertical distance between the beads.
Note that these refer to *center to center distances* (unlike the distance channel in Bluelake, which represents the bead surface to surface distance).
Note that all three parameters have to be specified in the same spatial unit (meters or micron).
The final parameter `is_y_oscillation` indicates whether the stage was oscillated in the y-direction.

The obtained correction factor can be used to correct the calibration factors::

Rd_corrected = factor * calibration["Rd"].value
Rf_corrected = calibration["Rf"].value / factor
stiffness_corrected = calibration["kappa"].value / factor**2

To correct a force trace, simply divide it by the correction factor::

corrected_force1x = f.force1x / factor

.. note::

This coupling model neglects effects from the surface. It is intended for measurements performed at the center of the flowcell.

.. note::

The model implemented here only supports beads that are aligned in the same plane.
It does not take a mismatch in the `z`-position of the beads into account.
In reality, the position in the focus depends on the bead radius and may be different for the two beads if they slightly differ in size :cite:`alinezhad2018enhancement` (Fig. 3).
At short bead-to-bead distances, such a mismatch would make the coupling less pronounced than the model predicts.

Fast Sensors
------------

Expand Down
50 changes: 50 additions & 0 deletions docs/tutorial/force_calibration/active_calibration.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
.. _bead_bead_tutorial:

Active calibration with two beads far away from the surface
-----------------------------------------------------------

.. warning::

The implementation of the coupling correction models is still alpha functionality.
While usable, this has not yet been tested in a large number of different scenarios.
The API can still be subject to change *without any prior deprecation notice*!
If you use this functionality keep a close eye on the changelog for any changes that may affect your analysis.

When performing active calibration with two beads, we get a lower fluid velocity around the beads than we would with a single bead.
This leads to a smaller voltage readout than expected and therefore a higher displacement sensitivity (microns per volt).
Failing to take this into account results in a bias.
Pylake offers a function to calculate a correction factor to account for the lower velocity around the bead.

.. note::

For more information on how these factors are derived, please refer to the :ref:`theory<bead_bead_theory>` section on this topic.

Appropriate correction factors for oscillation in x can be calculated as follows::

factor = lk.coupling_correction_2d(dx=5.0, dy=0, bead_diameter=bead_diameter, is_y_oscillation=False)

Here `dx` and `dy` represent the horizontal and vertical distance between the beads.
Note that these refer to *center to center distances* (unlike the distance channel in Bluelake, which represents the bead surface to surface distance).
Note that all three parameters have to be specified in the same spatial unit (meters or micron).
The final parameter `is_y_oscillation` indicates whether the stage was oscillated in the y-direction.

The obtained correction factor can be used to correct the calibration factors::

Rd_corrected = factor * calibration["Rd"].value
Rf_corrected = calibration["Rf"].value / factor
stiffness_corrected = calibration["kappa"].value / factor**2

To correct a force trace, simply divide it by the correction factor::

corrected_force1x = f.force1x / factor

.. note::

This coupling model neglects effects from the surface. It is intended for measurements performed at the center of the flowcell.

.. note::

The model implemented here only supports beads that are aligned in the same plane.
It does not take a mismatch in the `z`-position of the beads into account.
In reality, the position in the focus depends on the bead radius and may be different for the two beads if they slightly differ in size :cite:`alinezhad2018enhancement` (Fig. 3).
At short bead-to-bead distances, such a mismatch would make the coupling less pronounced than the model predicts.
2 changes: 1 addition & 1 deletion docs/tutorial/force_calibration/calibration_items.rst
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ To do this, we divide our data by the force sensitivity that was active at the s

To reproduce a calibration performed with Bluelake, the easiest way to extract all the relevant parameters is to use :meth:`~lumicks.pylake.force_calibration.calibration_item.ForceCalibrationItem.calibration_params()`::

>>> calibration_params = calibration.calibration_params()
>>> calibration_params = old_calibration.calibration_params()
... calibration_params
{'num_points_per_block': 2000,
'sample_rate': 78125,
Expand Down
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
77 changes: 74 additions & 3 deletions docs/tutorial/force_calibration/force_calibration.rst
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,18 @@ The most important parameters are:

- The bead diameter (in microns).
- The viscosity of the medium, which strongly depends on :ref:`temperature<temperature_theory>` which must also be provided. This viscosity of water at a particular temperature can be found using :func:`~lumicks.pylake.viscosity_of_water` :cite:`huber2009new`.
- The distance to the surface (in case of a surface experiment).

To find the viscosity of water at a particular temperature, Pylake uses :func:`~lumicks.pylake.viscosity_of_water` which implements the model presented in :cite:`huber2009new`.
When omitted, this function will automatically be used to look up the viscosity of water for that particular temperature

.. image:: figures/temperature_dependence.png
:nbattach:

.. note::

Note that for experiments that use a different medium than water, the viscosity at the experimental temperature should explicitly be provided.

Hydrodynamically correct model
""""""""""""""""""""""""""""""

Expand All @@ -45,17 +53,80 @@ Pylake uses values for water and polystyrene for the sample and bead density res
Experiments near the surface
""""""""""""""""""""""""""""

So far, we have only considered experiments performed deep in the flow-cell.
In reality, proximity of the flowcell surface to the bead leads to an increase in the drag force on the bead.

When doing experiments near the surface, it is recommended to provide a `distance_to_surface`.
This distance should be the distance from the center of the bead to the surface of the flowcell.
Since it can be challenging to determine this distance, it is recommended to use active calibration
when calibrating near the surface, since this makes calibration far less sensitive to mis-specification
of the bead diameter and height.

.. image:: figures/surface_calibration_workflow.png

Active or passive?
""""""""""""""""""

In the case of active calibration, it is mandatory to provide a nanostage signal, as well as a
guess of the driving frequency.
Active calibration has a few benefits.
When performing passive calibration, we base our calculations on a theoretical drag coefficient which depends on parameters that are only known with limited precision:

- The diameter of the bead :math:`d` in microns.
- The dynamic viscosity :math:`\eta` in Pascal seconds (which depends strongly on the local temperature around the bead, which is typically poorly known).
- The distance to the surface :math:`h` in microns.

For this reason, :ref:`Active calibration<active_calibration_theory>` is highly recommended when the bead is close to the surface, as it is less sensitive to the distance to the surface and bead diameter.
It is also recommended when the viscosity of the medium or the bead diameter are poorly known.
In the case of active calibration, it is mandatory to provide a nanostage signal, as well as a guess of the driving frequency.
Let's compare passive and active calibration for a surface experiment::

f = lk.File("test_data/near_surface_active_calibration.h5")

# Decalibrate the force data
volts = f.force1x / f.force1x.calibration[0].force_sensitivity

# Grab the nanostage data
driving_data = f["Nanostage position"]["X"]

# Since we will be doing a few calibrations, let's store the parameters in a dictionary
shared_parameters = {
"force_voltage_data": volts.data,
"bead_diameter": bead_diameter,
"temperature": 25,
"sample_rate": volts.sample_rate,
"driving_data": driving_data.data,
"driving_frequency_guess": 37,
"hydrodynamically_correct": False, # We will be too close to the surface for this model
}

Next, unpack this dictionary using the unpacking operator `**`::

>>> fit = lk.calibrate_force(
... **shared_parameters, active_calibration=True, distance_to_surface=distance_to_surface
... )
>>> print(fit["kappa"].value)
0.11662183772410809

And compare this to the passive calibration result::

>>> fit = lk.calibrate_force(
... **shared_parameters, active_calibration=False, distance_to_surface=distance_to_surface
... )
>>> print(fit["kappa"].value)
0.11763849764570819

These values are quite close.
However, if we do not provide the height above the surface, we can see that the passive calibration result suffers much more than the active calibration result (as passive calibration fully relies on a drag coefficient calculated from the physical input parameters)::

>>> print(lk.calibrate_force(**shared_parameters, active_calibration=False)["kappa"].value)
>>> print(lk.calibrate_force(**shared_parameters, active_calibration=True)["kappa"].value)
0.08616565751377737
0.11662183772410809

.. note::

Note that the drag coefficient `gamma_ex` that Pylake returns always corresponds to the drag coefficient extrapolated back to its bulk value.
This ensures that drag coefficients can be compared and carried over between experiments performed at different heights.
The field `local_drag_coefficient` contains an estimate of the local drag coefficient (at the provided height).

Axial force
"""""""""""
Expand Down Expand Up @@ -146,7 +217,7 @@ Excluding noise floors can be done by providing a list of tuples to the `exclude
plt.tight_layout()
plt.title(f"Stiffness: {calibration2.stiffness:.3f}, Force sensi: {calibration2.displacement_sensitivity:.2f}")

.. image:: figures/excluded_ranges.png
.. image:: figures/frequency_exclusion_ranges.png

Note that when plotting the calibration, we have used `show_excluded=True`, which shows the excluded ranges in the plot.
We can request these excluded ranges from the calibration item itself::
Expand Down
1 change: 1 addition & 0 deletions docs/tutorial/force_calibration/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -9,5 +9,6 @@ Find out about force calibration

calibration_items
force_calibration
active_calibration
low_level_api
robust_fitting
Loading

0 comments on commit 7a26ab6

Please sign in to comment.