Membrane thickness — analysis and visualization

This notebook covers post-pipeline analysis of the thickness and intensity profile files produced by cryocat.analysis.memthick (see tutorial Membrane thickness: run) via the memthick_analyze_plot module.

It assumes that the *_thickness.csv + *_int_profiles.pkl files exist in your outputs folder.

---

Table of contents

1. Pipeline

2. Setup

3. Boundary mode summary

4. Thickness distributions

5. 3D spatial visualization

6. Intensity profile summary

7. Binned intensity profiles

8. Motl export

9. Intensity-profile asymmetry

---

1 Pipeline

1.1 The three steps

The cryocat.analysis.memthick pipeline runs in three sequential steps. The first two are geometric; the third identifies minima and inflection points to calculateunderl membrane thickness from the udnerlying tomogram intensity profiles.

Step	Input	Output	What it computes
1 — Surface extraction	Segmentation MRC	*_vertices_normals.csv	Marching-cubes mesh of the membrane; surface normals; bilateral split into surface 1 / surface 2
2 — Geometric matching	Vertex CSV	*_matched_points.csv	One-to-one normal-cone pairing between the two surfaces; records the Euclidean distance along the normal (matched_points_distance_nm)
3 — Intensity profile analysis	Matched pairs + tomogram MRC	*_thickness.csv, *_int_profiles.pkl, *_boundary_stats.txt	Short 1-D intensity profiles extracted along each matched pair; detection of leaflet minima and outward maxima; inflection-point-based membrane thickness

Steps 1–2 are purely geometric — no tomogram is needed. Step 3 is where the membrane thickness is measured from the data.

---

1.2 What is membrane thickness here?

Two different distances appear in the output:

``matched_points_distance_nm`` (geometric)

The Euclidean distance between the matched segmentation surface points on surface 1 and surface 2, measured along the normal direction (hence, it’s a property of the segmentation).

``membrane_thickness_nm`` (profile-derived)

The distance between the inflection points of the intensity profile on the left and right sides of the membrane. An inflection point lies on the steepest part of the intensity ramp between a leaflet minimum and its outward maximum — it marks where the density gradient is largest and the points where we choose to define thickness. The distances between minima are also reported in case that’s the membrane thickness definition you’d prefer.

---

1.3 Intensity profile features

Each extracted profile runs along the vector connecting the matched point pair, centred at the geometric midpoint, and extends profile_half_width_nm in each direction.

The pipeline identifies five points per profile (where detectable):

Point	Symbol	Description
Left leaflet minimum	minima1	Lowest intensity on the left side (phosphate headgroups)
Right leaflet minimum	minima2	Same on the right side
Central maximum	central_max	Peak between the two minima; the hydrophobic core
Left outward maximum	—	Peak outside (left of) minima1 as the intensity returns to baseline
Right outward maximum	—	Peak outside (right of) minima2 as the intensity returns to baseline

The inflection point on each side lies on the slope between the leaflet minimum and the outward maximum — where intensity is changing fastest. Its position is found by a weighted-median slope analysis and refined by quadratic interpolation.

---

1.4 Detection modes — how detection_mode is assigned

Not every profile has a clear outward maximum on both sides. The pipeline handles this in two passes and records how each boundary was determined in the detection_mode column.

max_max — both boundaries from resolved outward maxima

Both the left and right outward maxima are found as genuine find_peaks peaks within anchor_search_nm of each leaflet minimum. The inflection point on each side is computed from the slope between the minimum and its maximum.

This is the highest-confidence assignment. Both sides of the membrane are well-resolved in the tomogram and the profile has a clear double-peak structure.

max_anchor / anchor_max — one boundary inferred by mirroring

One side has a clear outward maximum (max); the other does not. The missing side’s anchor position is predicted by mirroring the max-to-min span from the good side. The inflection point is then computed normally between the minimum and the mirrored anchor.

• max_anchor: left side has a real maximum; right side uses a mirrored anchor.

• anchor_max: right side has a real maximum; left side uses a mirrored anchor.

Both are moderate-confidence assignments. The membrane boundary is physically present but the intensity peak on one side is too weak or poorly separated to be detected as a strict maximum. The mirrored position uses symmetry to estimate where it should be.

In the analysis module, max_anchor and anchor_max are treated as equivalent and counted together. The distinction is preserved in the raw column for debugging.

minima_only — no inflection points; distance between minima

Neither outward maximum could be found. The pipeline falls back to reporting the distance between the two leaflet minima as minima_separation_nm.

• membrane_thickness_nm is set to NaN for these rows — they are different from inflection-point-based measurements.

• minima_separation_nm is always populated (for all rows, not just minima_only) and represents the inner minima-to-minima distance.

• minima_only rows appear in the exported CSV because they might still useful for spatial analysis.

Summary table:

detection_mode	Left boundary	Right boundary	membrane_thickness_nm	minima_separation_nm
max_max	Inflection point (real max)	Inflection point (real max)	✓	✓
max_anchor	Inflection point (real max)	Inflection point (mirrored)	✓	✓
anchor_max	Inflection point (mirrored)	Inflection point (real max)	✓	✓
minima_only	—	—	NaN	✓

---

1.5 Output files

File	Content
*_thickness.csv	One row per resolved profile; coordinates, membrane_thickness_nm (NaN for minima_only), minima_separation_nm, matched_points_distance_nm, detection_mode, surface normals
*_int_profiles.pkl	List of profile dicts (same order as rows in *_thickness.csv); each dict contains the raw intensity array, pixel size, matched-point geometry, and a features sub-dict with all detected landmarks
*_boundary_stats.txt	Text summary of boundary-finding statistics from the pipeline run

---

2 Setup

Point output_path and seg_base at your pipeline output folder.

membrane_names should match the label strings used when the pipeline was run — they are part of the output filenames.

%load_ext autoreload
%autoreload 2

from pathlib import Path
import memthick_analyze_plot as maple
from IPython.display import Image

output_path = Path("outputs")
seg_base    = "2140_z150to400_segmented"

membrane_names = ["ER", "IMM", "OMM"]

# Helper: build the path to a named output file for a given membrane
def out(membrane, suffix):
    return output_path / f"{seg_base}_{membrane}_{suffix}"

# Quick check that expected files exist
for membrane in membrane_names:
    for suffix in ["thickness.csv", "int_profiles.pkl"]:
        p = out(membrane, suffix)
        status = "OK" if p.exists() else "MISSING"
        print(f"{membrane:4s}  {suffix}  →  {status}")

ER    thickness.csv  →  OK
ER    int_profiles.pkl  →  OK
IMM   thickness.csv  →  OK
IMM   int_profiles.pkl  →  OK
OMM   thickness.csv  →  OK
OMM   int_profiles.pkl  →  OK

Load and analyse each membrane

load_membrane_data reads the thickness CSV and, when auto_discover_related_files=True, automatically finds the matching *_int_profiles.pkl and *_boundary_stats.txt in the same folder.

analyze_membrane_thickness computes summary statistics (mean, median, std, percentiles).

analyze_intensity_profiles summarises the extracted profiles and can optionally bin them by thickness.

# Pixel size
pixel_size_nm      = None     # None → auto-read from pickle

# Thickness filtering
outlier_removal    = None     # None |'percentile'

# Intensity profile axis
extension_range_nm = None     # None → inferred from profiles
                              # e.g. (-5.5, 5.5) to fix the x-axis window

membrane_data     = {}
thickness_results = {}
profile_results   = {}

for m in membrane_names:
    print(f"\n{'═'*60}")
    print(f"  {m}")
    print(f"{'═'*60}")

    # Load
    md = maple.load_membrane_data(
        thickness_csv=out(m, "thickness.csv"),
        auto_discover_related_files=True,
        pixel_size_nm=pixel_size_nm,
    )

    n_rows = len(md.thickness_df)
    n_prof = len(md.intensity_profiles)
    ps = (pixel_size_nm
          or (md.intensity_profiles[0].get("pixel_size") if md.intensity_profiles else None))
    ps_str = f"{ps:.4f} nm" if ps else "(unknown — pass pixel_size_nm)"
    print(f"  CSV rows : {n_rows:,}")
    print(f"  Profiles : {n_prof:,}")
    print(f"  Pixel size : {ps_str}")

    # Boundary mode counts
    info = maple.return_boundary_info(md)
    n_inf = info["n_finite_inflection_thickness_nm"]
    n_mo  = info["by_detection_mode"].get("minima_only", 0)
    print(f"\n  Inflection thickness (non-NaN): {n_inf:,}")
    print(f"  Minima-only rows (NaN thickness, min–min dist): {n_mo:,}")
    if info["by_detection_mode"]:
        stages = "  ".join(f"{k}: {v:,}" for k, v in info["by_detection_mode"].items())
        print(f"  detection_mode → {stages}")

    # Thickness analysis
    tr = maple.analyze_membrane_thickness(
        data=md,
        outlier_removal_method=outlier_removal,
    )
    s = tr.statistics
    n_fin = s.get("count_finite", s["count"])
    print(f"\n  Thickness (nm)  n_finite={n_fin:,}")
    print(f"    mean ± std : {s['mean']:.3f} ± {s['std']:.3f}")
    print(f"    median     : {s['median']:.3f}")
    print(f"    IQR        : {s['q25']:.3f} – {s['q75']:.3f}")
    print(f"    range      : {s['min']:.3f} – {s['max']:.3f}")

    # Intensity profile analysis
    pr = maple.analyze_intensity_profiles(
        data=md,
        extension_range_nm=extension_range_nm,
    )
    ext = pr.parameters.get("extension_range_nm", "?")
    pix = pr.parameters.get("pixel_size_nm")
    pix_str = f"{pix:.4f} nm" if pix else "?"
    print(f"\n  Profiles analysed : {len(pr.profiles):,}")
    if ext != "?":
        print(f"  Extension range   : {ext[0]:.2f} – {ext[1]:.2f} nm")
    print(f"  Pixel size used   : {pix_str}")

    membrane_data[m]     = md
    thickness_results[m] = tr
    profile_results[m]   = pr

════════════════════════════════════════════════════════════
  ER
════════════════════════════════════════════════════════════
Loaded 172110 intensity profiles from outputs/2140_z150to400_segmented_ER_int_profiles.pkl
  CSV rows : 172,110
  Profiles : 172,110
  Pixel size : 0.7840 nm

  Inflection thickness (non-NaN): 148,728
  Minima-only rows (NaN thickness, min–min dist): 23,382
  detection_mode → max_max: 81,261  max_anchor: 67,467  minima_only: 23,382

  Thickness (nm)  n_finite=148,728
    mean ± std : 5.781 ± 0.824
    median     : 5.671
    IQR        : 5.251 – 6.248
    range      : 2.391 – 8.000
Analyzing 172110 intensity profiles...

  Profiles analysed : 172,110
  Extension range   : -5.50 – 5.50 nm
  Pixel size used   : 0.7840 nm

════════════════════════════════════════════════════════════
  IMM
════════════════════════════════════════════════════════════
Loaded 678128 intensity profiles from outputs/2140_z150to400_segmented_IMM_int_profiles.pkl
  CSV rows : 678,128
  Profiles : 678,128
  Pixel size : 0.7840 nm

  Inflection thickness (non-NaN): 594,155
  Minima-only rows (NaN thickness, min–min dist): 83,973
  detection_mode → max_max: 371,655  max_anchor: 222,500  minima_only: 83,973

  Thickness (nm)  n_finite=594,155
    mean ± std : 5.544 ± 0.808
    median     : 5.371
    IQR        : 5.038 – 5.905
    range      : 2.095 – 8.000
Analyzing 678128 intensity profiles...

  Profiles analysed : 678,128
  Extension range   : -5.50 – 5.50 nm
  Pixel size used   : 0.7840 nm

════════════════════════════════════════════════════════════
  OMM
════════════════════════════════════════════════════════════
Loaded 495947 intensity profiles from outputs/2140_z150to400_segmented_OMM_int_profiles.pkl
  CSV rows : 495,947
  Profiles : 495,947
  Pixel size : 0.7840 nm

  Inflection thickness (non-NaN): 433,282
  Minima-only rows (NaN thickness, min–min dist): 62,665
  detection_mode → max_max: 260,979  max_anchor: 172,303  minima_only: 62,665

  Thickness (nm)  n_finite=433,282
    mean ± std : 5.530 ± 0.869
    median     : 5.278
    IQR        : 4.950 – 5.946
    range      : 2.200 – 8.000
Analyzing 495947 intensity profiles...

  Profiles analysed : 495,947
  Extension range   : -5.50 – 5.50 nm
  Pixel size used   : 0.7840 nm

---

3 Boundary mode summary

return_boundary_info reads the detection_mode column of the loaded CSV and reports:

• Total rows — all rows in the exported CSV (all resolved by construction).

• Inflection-point thickness rows — rows with a finite membrane_thickness_nm (max_max + max_anchor + anchor_max).

• Minima-only rows — rows where membrane_thickness_nm is NaN; only minima_separation_nm is available.

• Detection mode breakdown — count per mode.

for membrane in membrane_names:
    print(f"\n{'='*60}")
    print(f"  {membrane}")
    print(f"{'='*60}")
    info = maple.return_boundary_info(membrane_data[membrane])
    print("\n".join(info["lines"]))

============================================================
  ER
============================================================
Source: thickness_csv  |  Rows: 172,110
Inflection (non-NaN thickness): 148,728  (86.4%)
Minima-only: 23,382  (13.6%)
Detection mode:  max_max 81,261  |  max_anchor 67,467  |  minima_only 23,382

============================================================
  IMM
============================================================
Source: thickness_csv  |  Rows: 678,128
Inflection (non-NaN thickness): 594,155  (87.6%)
Minima-only: 83,973  (12.4%)
Detection mode:  max_max 371,655  |  max_anchor 222,500  |  minima_only 83,973

============================================================
  OMM
============================================================
Source: thickness_csv  |  Rows: 495,947
Inflection (non-NaN thickness): 433,282  (87.4%)
Minima-only: 62,665  (12.6%)
Detection mode:  max_max 260,979  |  max_anchor 172,303  |  minima_only 62,665

---

4 Thickness distributions

4.1 — Single membrane

plot_thickness_distribution shows a histogram (and optional KDE) of membrane_thickness_nm for one or more membranes.

Filtering parameters:

Parameter	Type	Description
thickness_range_nm	(min, max) or None	Display only rows with membrane_thickness_nm in this range. Rows with NaN thickness (minima_only) are dropped.
minima_separation_range_nm	(min, max) or None	Display only rows with minima_separation_nm in this range. Applied after thickness_range_nm (logical AND).
thickness_regime	str or list[str] or None	Keep only rows matching a detection_mode group: "max_max", "max_anchor", "minima_only", or a list like ["max_max", "max_anchor"].

Other parameters:

Parameter	Type	Description
histogram_bins	int	Number of histogram bins.
density_normalization	bool	Normalise y-axis to probability density (area = 1). Useful when comparing membranes with very different row counts.
show_mean_lines	bool	Draw vertical lines at the mean and median.
opacity	float	Bar opacity (0–1); useful when overlaying multiple membranes.

thickness_range_nm=(3.0, 8.0) #nm

fig = maple.plot_thickness_distribution(
    thickness_results["IMM"],
    membrane_names=["IMM"],
    thickness_range_nm=thickness_range_nm,
    histogram_bins=60,
    show_statistics=True,
    show_mean_lines=True,
    density_normalization=True,
)
#fig.show()
fig.write_image("plots/IMM_thickness_distribution.png", width=800, height=600)

Statistics for IMM:
  Mean:   5.54
  Std:    0.81
  Median: 5.37
  Count:  594,049

Image("plots/IMM_thickness_distribution.png")

4.2 Comparison of mutliple membranes

Pass a list of ThicknessAnalysisResults objects and a matching list of membrane_names to overlay distributions on the same axes.

fig = maple.plot_thickness_distribution(
    [thickness_results[m] for m in membrane_names],
    membrane_names=membrane_names,
    thickness_range_nm=thickness_range_nm,
    histogram_bins=60,
    density_normalization=True,
    show_mean_lines=True,
    opacity=0.6,
    plot_title="Membrane thickness comparison",
)
#fig.show()
fig.write_image("plots/thickness_distribution_comparison.png", width=800, height=600)

Statistics for ER:
  Mean:   5.78
  Std:    0.82
  Median: 5.67
  Count:  148,669

Statistics for IMM:
  Mean:   5.54
  Std:    0.81
  Median: 5.37
  Count:  594,049

Statistics for OMM:
  Mean:   5.53
  Std:    0.87
  Median: 5.28
  Count:  433,192

Image("plots/thickness_distribution_comparison.png")

4.3 — Leaflet min-to-min separation

plot_min_to_min_distribution plots the distance between the two detected minima (including minima_only rows where no other features were detected). You can disable this by setting thickness_regime.

fig = maple.plot_min_to_min_distribution(
    [thickness_results[m] for m in membrane_names],
    membrane_names=membrane_names,
    minima_separation_range_nm=(2.0, 7.0),
    histogram_bins=25,
    density_normalization=True,
    show_mean_lines=True,
    opacity=0.6,
    plot_title="Minima-to-minima distances",
)
#fig.show()
fig.write_image("plots/min_min_distances.png", width=800, height=600)

Statistics for ER:
  Mean:   4.76
  Std:    0.84
  Median: 4.68
  Count:  160,378

Statistics for IMM:
  Mean:   4.49
  Std:    0.86
  Median: 4.46
  Count:  642,327

Statistics for OMM:
  Mean:   4.49
  Std:    0.90
  Median: 4.46
  Count:  468,793

Image("plots/min_min_distances.png")

---

5 3D spatial visualization

These plots render the membrane in voxel space and colour each point by its local thickness or by surface identity.

5.1 — Thickness map in voxel space

plot_thickness_3d colour-codes each matched-point location by membrane_thickness_nm. Use sample_fraction to reduce the number of points for interactive performance (sample_fraction=0.3 keeps 30 % at random).

Filtering parameters:

Parameter	Type	Description
thickness_range_nm	(min, max) or None	Display only rows with membrane_thickness_nm in this range. Rows with NaN thickness (minima_only) are dropped.
minima_separation_range_nm	(min, max) or None	Display only rows with minima_separation_nm in this range. Applied after thickness_range_nm (logical AND).
thickness_regime	str or list[str] or None	Keep only rows matching a detection_mode group: "max_max", "max_anchor", "minima_only", or a list like ["max_max", "max_anchor"].

Other parameters:

Parameter	Type	Description
sample_fraction	float	Fraction of rows to display (random subsample).
marker_size	int	Point size in pixels.
color_scale	str	Any Plotly colorscale name ("OrRd", "Viridis", "RdBu", etc.).

fig = maple.plot_thickness_3d(
    [thickness_results[m] for m in membrane_names],
    membrane_names=membrane_names,
    thickness_range_nm=thickness_range_nm,
    sample_fraction=0.2, # helps with rendering
    marker_size=2,
    color_scale="OrRd",
)
#fig.show()
fig.write_image("plots/comparison_thickness_3d.png", width=1200, height=800)

Image("plots/comparison_thickness_3d.png")

5.2 — Surface 1 vs Surface 2 (point scatter or mesh)

plot_surfaces shows the two matched leaflet surfaces as coloured point clouds. By default, surface 1 (inner leaflet) is blue and surface 2 (outer leaflet) is orange.

When save_split_surface_meshes=True was passed to the pipeline, triangulated .ply meshes are written alongside the CSV. Pass ply_base_path to overlay the mesh.

fig = maple.plot_surfaces(
    thickness_results["OMM"],
    membrane_names=["OMM"],
    sample_fraction=0.3,
    marker_size=2,
)
#fig.show()
fig.write_image("plots/OMM_surface_scatter.png", width=1200, height=800)

OMM: Surface 1 points: 148,784, Surface 2 points: 148,784

Image("plots/OMM_surface_scatter.png")

fig = maple.plot_surfaces(
    thickness_results["OMM"],
    membrane_names=["OMM"],
    ply_base_path=output_path / seg_base,  # → outputs/2140_z150to400_segmented_OMM_surface1.ply
    mesh_opacity=0.5,
    show_scatter=False,  # mesh-only view; set True to see both
)
#fig.show()
fig.write_image("plots/OMM_surface_mesh.png", width=1200, height=800)

OMM: Surface 1 points: 495,947, Surface 2 points: 495,947

Image("plots/OMM_surface_mesh.png")

---

6 Intensity profile summary

plot_intensity_profile_summary stacks all (or a filtered subset of) profiles, interpolates them onto a common distance axis, and shows the mean profile with a 25th–75th percentile intensity band.

Vertical markers show where key points fall on average across all included profiles.

Marker display parameters:

Parameter	Default	What is shown
show_segmentation_boundary_markers	True	Where the geometrically matched segmentation surface points project onto the profile.
show_segmentation_boundary_distributions	True	IQR shadow (25th–75th percentile) around the positions.
show_inflection_point_markers	True	Mean left and right inflection point positions. Their separation equals the mean membrane_thickness_nm.
show_inflection_point_distributions	True	IQR shadow around the inflection point positions.
show_minima_midpoint	True	Geometric midpoint between minima.
show_minima	False	Mean positions of the left and right minima.
show_outward_maxima	False	Mean positions of the outward maxima (one per side). Useful for debugging boundary detection.

Filtering parameters:

Parameter	Type	Description
thickness_range_nm	(min, max) or None	Display only rows with membrane_thickness_nm in this range. Rows with NaN thickness (minima_only) are dropped.
minima_separation_range_nm	(min, max) or None	Display only rows with minima_separation_nm in this range. Applied after thickness_range_nm (logical AND).
thickness_regime	str or list[str] or None	Keep only rows matching a detection_mode group: "max_max", "max_anchor", "minima_only", or a list like ["max_max", "max_anchor"].

6.1 — Show all profiles

fig = maple.plot_intensity_profile_summary(
    profile_results["IMM"],
    membrane_names=["IMM"],
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=True,
    show_inflection_point_markers=True,
    show_inflection_point_distributions=True,
    show_minima_midpoint=True,
    plot_title="IMM intensity profiles summary",
)
#fig.show()
fig.write_image("plots/IMM_profiles_summary.png", width=800, height=600)

Image("plots/IMM_profiles_summary.png")

6.2 — Filter by thickness regime

Plotting profiles split by detection_mode is useful for quality control: if max_max and max_anchor profiles might look different (different mean thickness, different profile shape)

The thickness_regime parameter accepts any of the detection mode names or a list of them:

# Only profiles where both boundaries were detected from real maxima
thickness_regime="max_max"

# Only profiles where one boundary was inferred by mirroring (both left-real/right-mirrored and vice versa)
thickness_regime="max_anchor"

fig = maple.plot_intensity_profile_summary(
    profile_results["IMM"],
    membrane_names=["IMM"],
    thickness_regime="max_max",
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=True,
    show_inflection_point_markers=True,
    show_inflection_point_distributions=True,
    show_minima_midpoint=True,
    plot_title="IMM profiles with two resolved outward maxima",
)
#fig.show()
fig.write_image("plots/IMM_profiles_max_max.png", width=800, height=600)

Image("plots/IMM_profiles_max_max.png")

fig = maple.plot_intensity_profile_summary(
    profile_results["IMM"],
    membrane_names=["IMM"],
    thickness_regime="max_anchor",
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=True,
    show_inflection_point_markers=True,
    show_inflection_point_distributions=True,
    show_minima_midpoint=True,
    plot_title="IMM profiles with one resolved outward maxima and one inferred boundary",
)
#fig.show()
fig.write_image("plots/IMM_profiles_max_anchor.png", width=800, height=600)

Image("plots/IMM_profiles_max_anchor.png")

6.3 — Filter by thickness range

Combine thickness_range_nm with thickness_regime to isolate a specific population — for example, profiles in the 4–5 nm range.

fig = maple.plot_intensity_profile_summary(
    profile_results["IMM"],
    membrane_names=["IMM"],
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=True,
    show_inflection_point_markers=True,
    show_inflection_point_distributions=True,
    show_outward_maxima=True,
    show_minima=True,
    show_minima_midpoint=True,
    thickness_range_nm=(4.0, 5.0),
    plot_title="IMM profiles in the 4-5 nm thickness range",
)
#fig.show()
fig.write_image("plots/IMM_profiles_4_5nm.png", width=800, height=600)

Image("plots/IMM_profiles_4_5nm.png")

---

7 Binned intensity profiles

plot_intensity_profile_binned splits profiles into thickness bins and plots the mean profile per bin. This reveals whether the profile shape changes with membrane thickness — for example, thicker membranes might show a broader central maximum or more separated minima.

Binning methods:

binning_method	Description
"quantile" (default)	Equal number of profiles per bin. Each bin has the same statistical weight.
"equal_width"	Equal nm range per bin. Bins near the mode of the distribution will contain far more profiles than tails.
Custom list	Pass [(min1, max1, "label1"), (min2, max2, "label2"), ...] to define arbitrary windows.

Key parameters:

Parameter	Type	Description
thickness_bins	int or list or None	Number of auto-bins, explicit bin list, or None (defaults to 4 quartiles).
binning_method	str	"quantile" or "equal_width". Ignored when a custom bin list is passed.

7.1 — Equal number of profiles per bin

fig = maple.plot_intensity_profile_binned(
    profile_results["IMM"],
    membrane_names=["IMM"],
    thickness_bins=4,
    binning_method="quantile",
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=False, # too crowded
    show_inflection_point_markers=True,
    show_inflection_point_distributions=False, # too crowded
    show_minima_midpoint=False,
    plot_title="IMM profiles binned by thickness quartiles",
)
#fig.show()
fig.write_image("plots/IMM_profiles_binned_quartiles.png", width=800, height=600)

Using 4 bins (thickness_nm) for IMM:
  Bin 1 (infl.): 2.1-5.0 nm: 148,539 profiles
  Bin 2 (infl.): 5.0-5.4 nm: 148,539 profiles
  Bin 3 (infl.): 5.4-5.9 nm: 148,539 profiles
  Bin 4 (infl.): 5.9-8.0 nm: 148,539 profiles

Image("plots/IMM_profiles_binned_quartiles.png")

7.2 — Equal-width bins

fig = maple.plot_intensity_profile_binned(
    profile_results["IMM"],
    membrane_names=["IMM"],
    thickness_bins=4,
    binning_method="equal_width",
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=False, # too crowded
    show_inflection_point_markers=True,
    show_inflection_point_distributions=False, # too crowded
    show_minima_midpoint=False,
    plot_title="IMM profiles binned by equal width thickness ranges",
)
#fig.show()
fig.write_image("plots/IMM_profiles_binned_equal_width.png", width=800, height=600)

Using 4 bins (thickness_nm) for IMM:
  Bin 1 (infl.): 2.1-3.6 nm: 1,256 profiles
  Bin 2 (infl.): 3.6-5.0 nm: 150,493 profiles
  Bin 3 (infl.): 5.0-6.5 nm: 367,124 profiles
  Bin 4 (infl.): 6.5-8.0 nm: 75,282 profiles

Image("plots/IMM_profiles_binned_equal_width.png")

7.3 — Custom thickness windows

Pass explicit (min_nm, max_nm, label) tuples to define your own bins.

fig = maple.plot_intensity_profile_binned(
    profile_results["IMM"],
    membrane_names=["IMM"],
    thickness_bins=[
        (0.0, 4.0,  "< 4 nm"),
        (4.0, 5.0,  "4 – 5 nm"),
        (5.0, 6.0,  "5 – 6 nm"),
        (6.0, 20.0, "> 6 nm"),
    ],
    show_segmentation_boundary_markers=True,
    show_segmentation_boundary_distributions=False, # too crowded
    show_inflection_point_markers=True,
    show_inflection_point_distributions=False, # too crowded
    show_minima_midpoint=False,
    plot_title="IMM profiles binned by custom thickness ranges",
)
#fig.show()
fig.write_image("plots/IMM_profiles_binned_custom.png", width=800, height=600)

Using 4 bins (thickness_nm) for IMM:
  < 4 nm: 5,753 profiles
  4 – 5 nm: 126,149 profiles
  5 – 6 nm: 325,153 profiles
  > 6 nm: 137,398 profiles

Image("plots/IMM_profiles_binned_custom.png")

---

8 Motl export

save_thickness_motls writes motive-list .em files (one for each surface) containing the XYZ coordinates, Euler angles derived from the surface normals, and a thickness value as the score column.

These files can be opened in ChimeraX for spatial visualization of the thickness map and overlaid with the tomogram.

Key parameters:

Parameter	Type	Description
thickness_csv	path	Path to the *_thickness.csv output from Stage 3.
score_column	str	Column written to the motl score field. Default "membrane_thickness_nm" exports inflection-point distances. Pass "minima_separation_nm" to export minima distances.
sample_fraction	float	Write a subsampled motl alongside the full one — useful for faster visualization in ChimeraX without loading hundreds of thousands of subtomograms.
output_path	path	Directory for the output .em files.

for m in membrane_names:
    maple.save_thickness_motls(
        thickness_csv=out(m, "thickness.csv"),
        sample_fraction=0.20, # for easier visualization in ChimeraX
        output_path=output_path,
        score_column="membrane_thickness_nm",
    )
    print(f"{m}: done")

Saved full motls with 148728 points
Saved subsampled motls with 29745 points
ER: done
Saved full motls with 594155 points
Saved subsampled motls with 118831 points
IMM: done
Saved full motls with 433282 points
Saved subsampled motls with 86656 points
OMM: done

---

9 Intensity-profile asymmetry

analyze_intensity_asymmetry quantifies whether the two leaflet minima intenisites systematically differ compared to the central maximum in each profile. This is a profile-shape metric, not a direct measurement of compositional leaflet asymmetry.

How the score is defined

All profiles within a thickness bin (or spatial neighbourhood) are first aggregated by taking the median intensity at the left minimum, right minimum, and central maximum. The asymmetry score is then computed once from these aggregated values to suppresses noise from poorly resolved individual profiles.

The two intensity differences from the aggregated profile are:

• Left side: intensity at the central maximum minus intensity at the left minimum

• Right side: intensity at the central maximum minus intensity at the right minimum

The unity-based asymmetry score is:

score = max(left_side, right_side) / min(left_side, right_side)

• Score 1.0 = the two minima have equal intensities compared to the central maximum (symmetric profile)

• Score > 1.0 = one side’s minima intensity differs; the ratio tells you by how much

• Percentage = (score − 1) × 100

Two binning modes

binning_mode	How profiles are grouped	What the score measures
"thickness"	Equal-width bins along membrane_thickness_nm	Does intensity-profile shape change with thickness?
"vicinity"	All profiles within a spatial radius around each seed voxel	Where on the membrane surface do asymmetric profiles cluster?

Bins or seed points with fewer than min_profiles_per_bin valid profiles are rejected.

Interpreting results — general guidance

High asymmetry scores can arise from:

1. Missing-wedge artefacts

2. Poorly resolved membranes

3. Genuine variation — in some cases, real compositional or curvature-related differences may produce asymmetric profiles; but this interpretation requires corroboration from independent evidence and exclusion of the two points above.

Always run ``plot_asymmetry_detection_examples`` to visually inspect the aggregated median profiles for a representative set of bins. If the aggregated profile looks malformed, the asymmetry score is unreliable.

Consider: - Re-running with a different neighbourhood radius (vicinity_radius_nm) or bin size (thickness_bin_size_nm) - Enabling outlier_removal_method="percentile" (with outlier_percentile_range=(5, 95) or narrower) to suppress bins driven by extreme outliers

9.1 — Thickness-based binning

Groups profiles into equal-width bins of thickness_bin_size_nm nm and computes one asymmetry score per bin.

Because profiles are grouped by a global thickness value rather than spatial proximity, this mode is relatively conservative: each bin pools a large, spatially distributed set of profiles, which averages out local artefacts and tends to produce fewer extreme scores. The trade-off is that spatial information is lost — a high score in a thickness bin tells you that profiles of that thickness are asymmetric on average, but not where on the membrane surface they are.

Key parameters:

Parameter	Type	Description
thickness_bin_size_nm	float	Width of each thickness bin in nm (default 0.5).
min_profiles_per_bin	int	Bins with fewer valid profiles are rejected.
thickness_range_nm	(min, max) or None	Restrict analysis to this thickness window.
outlier_removal_method	None \| "percentile"	Filter outlier bins from the results. "percentile" removes bins whose score falls outside outlier_percentile_range.
outlier_percentile_range	(float, float)	Lower and upper percentile cutoffs used when outlier_removal_method="percentile" (default (5, 95)).

# Asymmetry analysis parameters
asym_thickness_bin_size_nm = 0.1
asym_min_profiles_per_bin  = 50

asymmetry_results_thickness = {}

for m in membrane_names:
    print(f"\n{'═'*60}")
    print(f"  {m}  (thickness mode)")
    print(f"{'═'*60}")

    ar = maple.analyze_intensity_asymemtry(
        data=membrane_data[m],
        binning_mode="thickness",
        thickness_bin_size_nm=asym_thickness_bin_size_nm,
        min_profiles_per_bin=asym_min_profiles_per_bin,
        thickness_range_nm=thickness_range_nm,
        outlier_removal_method=outlier_removal,
    )

    s = ar.overall_statistics
    print(f"  Valid bins       : {s.get('valid_bins', '?')}")
    print(f"  Rejected bins    : {len(ar.rejected_bins)}")
    print(f"  Mean score       : {s.get('mean_asymmetry', float('nan')):.3f}")
    print(f"  Median score     : {s.get('median_asymmetry', float('nan')):.3f}")
    print(f"  Mean percentage  : {s.get('mean_asymmetry_percent', float('nan')):.1f} %")

    asymmetry_results_thickness[m] = ar

════════════════════════════════════════════════════════════
  ER  (thickness mode)
════════════════════════════════════════════════════════════
Starting asymmetry analysis on 172110 profiles (binning_mode='thickness')...
Filtered to thickness range 3.0-8.0 nm: 148669 profiles remain
Created 50 thickness bins from 3.00 to 8.00 nm
Successfully computed asymmetry for 46 bins/seeds
Rejected: 4
  Valid bins       : 46
  Rejected bins    : 4
  Mean score       : 1.033
  Median score     : 1.022
  Mean percentage  : 3.3 %

════════════════════════════════════════════════════════════
  IMM  (thickness mode)
════════════════════════════════════════════════════════════
Starting asymmetry analysis on 678128 profiles (binning_mode='thickness')...
Filtered to thickness range 3.0-8.0 nm: 594049 profiles remain
Created 50 thickness bins from 3.00 to 8.00 nm
Successfully computed asymmetry for 50 bins/seeds
Rejected: 0
  Valid bins       : 50
  Rejected bins    : 0
  Mean score       : 1.024
  Median score     : 1.013
  Mean percentage  : 2.4 %

════════════════════════════════════════════════════════════
  OMM  (thickness mode)
════════════════════════════════════════════════════════════
Starting asymmetry analysis on 495947 profiles (binning_mode='thickness')...
Filtered to thickness range 3.0-8.0 nm: 433192 profiles remain
Created 50 thickness bins from 3.00 to 8.00 nm
Successfully computed asymmetry for 49 bins/seeds
Rejected: 1
  Valid bins       : 49
  Rejected bins    : 1
  Mean score       : 1.049
  Median score     : 1.043
  Mean percentage  : 4.9 %

Distribution of asymmetry scores

fig = maple.plot_asymmetry_distribution(
    [asymmetry_results_thickness[m] for m in membrane_names],
    membrane_names=membrane_names,
    asymmetry_range = (1.0, 25.0),
    use_percentages=True,       # True → show (score−1)×100 %
    plot_title="Intensity profile asymmetry (binned by thickness)",
)
#fig.show()
fig.write_image("plots/asymmetry_distribution_thickness.png", width=800, height=600)

Image("plots/asymmetry_distribution_thickness.png")

3D spatial map of asymmetry scores

fig = maple.plot_asymmetry_3d(
    asymmetry_results_thickness["IMM"],
    membrane_names=["IMM"],
    asymmetry_range = (1.0, 25.0),
    sample_fraction=0.8,
    marker_size=4,
    color_scale="RdBu_r",
    plot_title="IMM intensity profile asymmetry (3D map, binned by thickness)",
)
#fig.show()
fig.write_image("plots/IMM_asymmetry_3d_thickness_binned.png", width=1200, height=800)

Image("plots/IMM_asymmetry_3d_thickness_binned.png")

Aggregated profile examples

plot_asymmetry_detection_examples shows the aggregated median profile for a selection of bins or seed points, with the two minima and central maximum marked. This is the most important quality-control step: inspect whether the marked features are consistent with a well-resolved double-peak profile.

Only bins for which profile data was cached during analysis are available (controlled by max_cached_profile_examples in analyze_intensity_asymmetry).

fig = maple.plot_asymmetry_detection_examples(
    asymmetry_results_thickness["IMM"],
    number_of_examples=9,          # number of bins to show; evenly sampled across score range
    plot_title="IMM — asymmetric bin examples",
)
#fig.show()
fig.write_image("plots/IMM_asymmetry_examples_thickness_binned.png", width=1200, height=800)

Image("plots/IMM_asymmetry_examples_thickness_binned.png")

9.2 — Neighborhood-based binning

Instead of grouping by thickness, profiles are grouped spatially: for each seed voxel on surface 1, all profiles within neighborhood_radius_nm are pooled and one asymmetry score is computed from their median. The result is a dense spatial map of intensity-profile asymmetry across the membrane surface.

Interpret neighborhood-mode results with caution.

Unlike thickness-based binning, each seed point produces its own score from a spatially local pool of profiles. This means artefacts that are spatially correlated (e.g. the missing wedge) will appear as coherent high-asymmetry regions in the map. A cluster of high scores does not necessarily indicate genuine asymmetry; it may simply mark a region where the tomogram is poorly sampled or reconstructed.

Before interpreting spatial patterns: - Display the 3D asymmetry map alongside the missing-wedge orientation of your tilt series and check whether high-asymmetry regions align with the expected missing-wedge shadow. - Run plot_asymmetry_detection_examples on a representative selection of high-scoring seed points and verify that the aggregated profiles look well-resolved and double-peaked. - Consider re-running with a different neighborhood_radius_nm. A larger radius pools more profiles per seed (more robust scores, coarser spatial resolution). - If there is a long tail of extreme scores, apply outlier_removal_method="percentile" to suppress the most extreme bins before visualisation.

Key parameters:

Parameter	Type	Description
neighborhood_radius_nm	float	Radius of the spatial neighbourhood in nm.
neighborhood_seed_fraction	float	Fraction of surface-1 voxels used as seeds (1.0 = all). Reduce for speed on large membranes.
min_profiles_per_bin	int	Seeds with fewer valid neighbours are skipped.
outlier_removal_method	None \| "percentile"	Filter outlier seed scores from the results.
outlier_percentile_range	(float, float)	Percentile cutoffs for outlier removal (default (5, 95)).
random_seed	int	Random seed for reproducible subsampling when neighborhood_seed_fraction < 1.0.

asymmetry_results_neighborhood = {}

for m in membrane_names:
    print(f"\n{'═'*60}")
    print(f"  {m}  (neighborhood mode)")
    print(f"{'═'*60}")

    ar = maple.analyze_intensity_asymemtry(
        data=membrane_data[m],
        binning_mode="neighborhood",
        neighborhood_radius_nm=20.0,
        neighborhood_seed_fraction=0.1,   # use 10 % of voxels as seeds for speed
        min_profiles_per_bin=asym_min_profiles_per_bin,
        thickness_range_nm=thickness_range_nm,
        outlier_removal_method=outlier_removal,
        random_seed=42,
    )

    s = ar.overall_statistics
    print(f"  Valid seed points : {s.get('valid_bins', '?')}")
    print(f"  Rejected seeds    : {len(ar.rejected_bins)}")
    print(f"  Mean score        : {s.get('mean_asymmetry', float('nan')):.3f}")
    print(f"  Median score      : {s.get('median_asymmetry', float('nan')):.3f}")
    print(f"  Mean percentage   : {s.get('mean_asymmetry_percent', float('nan')):.1f} %")

    asymmetry_results_neighborhood[m] = ar

════════════════════════════════════════════════════════════
  ER  (neighborhood mode)
════════════════════════════════════════════════════════════
Starting asymmetry analysis on 172110 profiles (binning_mode='neighborhood')...
Filtered to thickness range 3.0-8.0 nm: 148669 profiles remain
Building KDTree on 148669 surface-1 points, radius=20.00 nm (25.51 vox), 14866 seeds (seed_fraction=0.1)...
Querying neighborhoods (bulk)...
Caching median profiles for 200 representative seeds...
Successfully computed asymmetry for 14866 bins/seeds
Rejected: 0
  Valid seed points : 14866
  Rejected seeds    : 0
  Mean score        : 1.076
  Median score      : 1.062
  Mean percentage   : 7.6 %

════════════════════════════════════════════════════════════
  IMM  (neighborhood mode)
════════════════════════════════════════════════════════════
Starting asymmetry analysis on 678128 profiles (binning_mode='neighborhood')...
Filtered to thickness range 3.0-8.0 nm: 594049 profiles remain
Building KDTree on 594049 surface-1 points, radius=20.00 nm (25.51 vox), 59404 seeds (seed_fraction=0.1)...
Querying neighborhoods (bulk)...
Caching median profiles for 200 representative seeds...
Successfully computed asymmetry for 59404 bins/seeds
Rejected: 0
  Valid seed points : 59404
  Rejected seeds    : 0
  Mean score        : 1.074
  Median score      : 1.061
  Mean percentage   : 7.4 %

════════════════════════════════════════════════════════════
  OMM  (neighborhood mode)
════════════════════════════════════════════════════════════
Starting asymmetry analysis on 495947 profiles (binning_mode='neighborhood')...
Filtered to thickness range 3.0-8.0 nm: 433192 profiles remain
Building KDTree on 433192 surface-1 points, radius=20.00 nm (25.51 vox), 43319 seeds (seed_fraction=0.1)...
Querying neighborhoods (bulk)...
Caching median profiles for 200 representative seeds...
Successfully computed asymmetry for 43319 bins/seeds
Rejected: 0
  Valid seed points : 43319
  Rejected seeds    : 0
  Mean score        : 1.081
  Median score      : 1.063
  Mean percentage   : 8.1 %

Distribution of asymmetry scores (neighborhood mode)

fig = maple.plot_asymmetry_distribution(
    [asymmetry_results_neighborhood[m] for m in membrane_names],
    membrane_names=membrane_names,
    use_percentages=True,
    asymmetry_range = (1.0, 25.0),
    plot_title="Intensity profile asymmetry (binned by spatial proximity)",
)
#fig.show()
fig.write_image("plots/asymmetry_distribution_neighborhood.png", width=800, height=600)

Image("plots/asymmetry_distribution_neighborhood.png")

3D spatial map (neighborhood mode)

fig = maple.plot_asymmetry_3d(
    asymmetry_results_neighborhood["IMM"],
    membrane_names=["IMM"],
    sample_fraction=0.8,
    asymmetry_range = (1.0, 25.0),
    marker_size=4,
    color_scale="RdBu_r",
    plot_title="IMM intensity profile asymmetry (3D map, binned by spatial proximity)",
)
#fig.show()
fig.write_image("plots/IMM_asymmetry_3d_neighborhood.png", width=1200, height=800)

Image("plots/IMM_asymmetry_3d_neighborhood.png")

fig = maple.plot_asymmetry_detection_examples(
    asymmetry_results_neighborhood["IMM"],
    number_of_examples=9,          # number of bins to show; evenly sampled across score range
    plot_title="IMM — asymmetric bin examples (binned by spatial proximity)",
)
#fig.show()
fig.write_image("plots/IMM_asymmetry_examples_neighborhoold_binned.png", width=1200, height=800)

Image("plots/IMM_asymmetry_examples_neighborhoold_binned.png")