Anisotropic stresses#
This notebook demonstrates the calculation of anisotropic stresses of injected oil droplets. The measure properties included the following:
\(\sigma^A(\vec{x}_s)\): anisotropic stress on surface
\(\sigma^A_T(\vec{x}_s)\): anisotropic stress on fitted ellipsoid (e.g. tissue-level)
import napari
import napari_stress
from napari_stress import measurements, approximation
import matplotlib.pyplot as plt
import numpy as np
napari.manifest -> 'napari' could not be imported: Could not find file 'builtins.yaml' in module 'napari'
config#
We need to set some parameters:
max_degree
: Degree of the spherical harmonics expansion. See also glossary.\(\gamma\)-value: Interfacial surface tension of an injected oil droplet in \(\frac{mN}{m}\)
max_degree = 5
gamma = 26.0
# Get some sample data
pointcloud = napari_stress.get_droplet_point_cloud()[0]
viewer = napari.Viewer(ndisplay=3)
napari.manifest -> 'napari' could not be imported: Could not find file 'builtins.yaml' in module 'napari'
layer_raw = viewer.add_layer(napari.layers.Layer.create(data = pointcloud[0], meta=pointcloud[1], layer_type=pointcloud[2]))
napari.utils.nbscreenshot(viewer)
Ellipsoid fit#
First, we fit an ellipsoid to the pointcloud. The ellipsoid represents the tissue-scale stresses.
expander = approximation.EllipsoidExpander()
expander.fit(pointcloud[0][:, 1:])
fitted_ellipsoid = expander.coefficients_
fitted_ellipse_points = expander.fit_expand(pointcloud[0][:, 1:])
We then measure the curvatures on every point on the surface:
curvature_on_ellipsoid = measurements.curvature_on_ellipsoid(fitted_ellipsoid, fitted_ellipse_points)
H_ellipsoid_major_medial_minor = measurements.mean_curvature_on_ellipse_cardinal_points(fitted_ellipsoid)
H_ellipsoid_major_medial_minor
[0.08211623521117135, 0.06563346684868424, 0.055253239527228264]
layer_raw.visible = False
ellipsoid_layer = viewer.add_vectors(fitted_ellipsoid)
viewer.add_points(fitted_ellipse_points, features = curvature_on_ellipsoid[1]['features'], face_color = 'mean_curvature', size=0.5, face_colormap = 'jet')
napari.utils.nbscreenshot(viewer)
We now do a spherical harmonics expansion of the ellipsoid points. From this, we can then calculate the average mean curvature \(H_0\) and the mean curvature \(H_i\) at every point on the surface.
fitted_ellipse_points_spherical_harmoniocs = napari_stress.fit_spherical_harmonics(fitted_ellipse_points, max_degree=max_degree)
coefficients = fitted_ellipse_points_spherical_harmoniocs[1]['metadata']['spherical_harmonics_coefficients']
quadrature_points, lebedev_fit = napari_stress.lebedev_quadrature(coefficients=coefficients, number_of_quadrature_points=500, use_minimal_point_set=False)
manifold_ellipsoid = napari_stress.create_manifold(quadrature_points, lebedev_fit=lebedev_fit, max_degree=max_degree)
Note: \(H_0\) is calculated using the spherical-harmonics related curvature definition. The point-wise curvature \(H_i\) is derived by projecting the quadrature points on the surface of the ellipsoid and using the respective measurement function:
_, _, H0_ellipsoid = measurements.calculate_mean_curvature_on_manifold(manifold_ellipsoid)
quadrature_points_on_ellipsoid = approximation.expand_points_on_ellipse(fitted_ellipsoid, quadrature_points)
mean_curvature_ellipsoid = measurements.curvature_on_ellipsoid(fitted_ellipsoid, quadrature_points_on_ellipsoid)[1]['features']['mean_curvature']
Spherical harmonics expansion#
We now do a spherical harmonics expansion of the raw pointcloud. This also captures the higher orders of curvature on the surface of the droplet.
fitted_points = napari_stress.fit_spherical_harmonics(pointcloud[0], max_degree=max_degree)
fitted_points_layer = viewer.add_points(fitted_points[0], **fitted_points[1], name='spherical harmonics expansion')
coefficients = fitted_points[1]['metadata']['spherical_harmonics_coefficients']
quadrature_points, lebedev_fit = napari_stress.lebedev_quadrature(coefficients=coefficients, number_of_quadrature_points=500, use_minimal_point_set=False)
manifold_droplet = napari_stress.create_manifold(quadrature_points, lebedev_fit=lebedev_fit, max_degree=max_degree)
mean_curvature_spherical_harmonics, _, H0_surface_integral_spherical_harmonics = measurements.calculate_mean_curvature_on_manifold(manifold_droplet)
Stresses#
We first calculate the tissue stress tensors along the ajor axes of the ellipsoid and the cartesian cardinal directions (x/y/z):
ellptical, cartesian = measurements.tissue_stress_tensor(fitted_ellipsoid, H0_ellipsoid, gamma=gamma)
Secondly, we calculate the anisotropic stress for every point on the surface:
stress, stress_tissue, stress_droplet = measurements.anisotropic_stress(mean_curvature_spherical_harmonics, H0_surface_integral_spherical_harmonics,
mean_curvature_ellipsoid, H0_ellipsoid,
gamma=gamma)
nbins = 20
fig, axes = plt.subplots(ncols=3, figsize=(15,5))
axes[0].hist(stress, nbins)
axes[0].set_title('Stress')
axes[1].hist(stress_tissue, nbins)
axes[1].set_title('tissue-scale stresses')
axes[2].hist(stress_droplet, nbins)
axes[2].set_title('Cell scale stresses')
for ax in axes:
ax.set_xlabel('anisotropic stress [mN/m]')
ax.set_ylabel('Occurrences [#]')
Lastly, we can also calculate the anisotropic tissue stressses \(\sigma^A_T\) from the fitted ellipsoid above:
ansiostropic_tissue_stress = measurements.maximal_tissue_anisotropy(fitted_ellipsoid, gamma=gamma)
ansiostropic_tissue_stress
1.396875775565042