#!/usr/bin/env python3 """Tier-2 item 16 -- anisotropic osmotic shrinkage and the axis-ratio prediction. Why this run matters -------------------- The model imposes ISOTROPIC natural shrinkage: `_target_axis_scales` returns `volume_ratio ** (1/3)` on all three axes unless `axis_scale_factors` is set, and `_growth_tensor` then raises that vector to the load factor (to the load factor times `shell_eigenstrain_fraction` in the shell). Adherent-cell nuclei under hyperosmotic stress do not do that: they flatten preferentially in height rather than contracting isotropically, and anisotropic growth redistributes the shell resultant instead of scaling it. Whether the reported unloading survives a realistic flattening mode has never been tested; the epistemic UQ varied the natural volume ratio but pinned the shape (`axis_scale_factors` is set to `[isotropic_scale] * 3` in run_reynolds_finan_epistemic_uq.py). The second reason is cheaper and better. The predicted change in the three semi-axis ratios is a NON-TRIVIAL directly measurable output. The volume ratio is an INPUT to this model -- it is prescribed from osmolarity and solid fraction -- so agreement on volume is not evidence. The shape change is not prescribed: it comes out of the shell/core mechanics. Confocal z-stacks before and after a 580 mOsm step cost a fraction of the proposed Flipper-TR FLIM experiment and measure exactly this. The axis-ratio prediction currently appears in one sentence of the wet-lab section, in no figure, and in no row of the claim matrix. What is imposed (state this verbatim in the manuscript) ------------------------------------------------------- One scalar `--anisotropy` a. With natural volume ratio V (V < 1 for shrinkage), the per-axis natural stretches are a partition of the SAME volume ratio in the exponent: s_x = V ** ((1 - a) / 3) s_y = V ** ((1 - a) / 3) s_z = V ** ((1 + 2 a) / 3) so that, exactly and for every a, s_x * s_y * s_z = V ** ((1-a)/3 + (1-a)/3 + (1+2a)/3) = V, and the imposed natural height-to-width stretch ratio is s_z / s_x = V ** a. a = 0 is the published isotropic case. a = 1 puts the entire natural volume change into the short z axis at fixed natural footprint. a > 1 would require in-plane natural EXPANSION and is rejected unless --allow-extreme is passed. Because the Control configuration has V = 1, every s_i is exactly 1 there and the Control problem is untouched at any a: the sweep perturbs the hyperosmotic natural state only. The shell inherits the same anisotropy direction with exponent `shell_eigenstrain_fraction`, which is the existing coupling, not a new assumption. The summary CSV column `natural_volume_ratio` is `prod(target_axis_scales)`, so it must stay at its published value across the whole sweep. The script checks that and refuses to report a row where it moved. How the axis ratios are obtained -------------------------------- The run outputs do NOT expose the deformed semi-axes: `_deformed_geometry_summary` records outer volume and outer area only, and `_summary_rows` carries `realized_outer_volume_ratio` and `realized_outer_area_ratio` and nothing about shape. Volume and area alone cannot separate a flattened ellipsoid from an isotropically shrunken one, so the CSV is insufficient by construction. They can still be recovered without touching upstream code, because each case directory stores the vertex displacement field: /nonlinear_nodal_displacement.npy (n_points, 3) /case_config.json (records what was imposed) Deformed coordinates are `mesh.points + displacement`, and this script reports two independent estimates of the deformed semi-axes: * half-extent a_i = (max - min) / 2 along each axis. The outer surface bounds the body, so the extremes are outer-surface points. * inertia-equivalent ellipsoid, from the exact per-tetrahedron second moments of the whole deformed body: for a solid ellipsoid the volume covariance is diag(a^2, b^2, c^2) / 5, so a_i = sqrt(5 * eigenvalue). This uses every element rather than four extreme vertices and is the number to quote. Both need the mesh, so they need meshio. For P2 runs the stored displacement is nodal at the mesh vertices, so the geometry is evaluated as if linear; the semi-axes are accurate to the vertex geometry, not to the P2 surface. The mesh is a polyhedron inscribed in the ellipsoid, so the equivalent semi-axes are biased low by O(h^2) -- about 0.3 % on a 1280-facet test sphere, and less on these meshes. The bias is identical in both conditions on the same mesh and cancels in the Control-to-hyperosmotic ratio, which is the reason the ratio CHANGE is the reported quantity and the absolute semi-axes are context. Required code change (only to put the numbers in the source data) ----------------------------------------------------------------- This script needs no upstream edit. Publishing the prediction does: 1. nonlinear_solver.py, `_deformed_geometry_summary` (and the same block at the end of `_write_outer_surface_resultants`): add `reference_outer_semi_axes_um` and `deformed_outer_semi_axes_um`, computed from `points` and `points + displacement` over the outer facets exactly as `_closed_surface_volume_area` already walks them. 2. run_reynolds_finan_heterogeneous_two_group.py, `_summary_rows`: add `deformed_outer_semi_axis_x_um`, `_y_um`, `_z_um` from that summary. Three floats per case, no solver change, and the measurable prediction becomes part of the shipped source data instead of a post-hoc reconstruction. How the anisotropy is injected ------------------------------ `run_analysis` builds its case dicts through `build_template_case_configs`, and those are plain JSON-shaped dicts: `_add_heterogeneous_provenance` calls `config.setdefault(...)` on them. The upstream epistemic-UQ runner already writes `hyper["stimulus"].update({"volume_ratio": ..., "axis_scale_factors": [...]})`, so the key is part of the supported schema. This script therefore wraps `runner.build_template_case_configs` -- the same runtime-patch discipline `_common.patch_options` uses for `_options` -- and leaves the upstream tree untouched. If that hook ever stops returning dicts, the script stops and prints the exact change needed rather than faking a result. Usage: python3 16_anisotropic_growth_axis_ratio.py --project-root /path/to/tree \ --out-root /path/to/out/anisotropy --anisotropy 0 0.25 0.5 0.75 1.0 --dry-run """ from __future__ import annotations import argparse import importlib import json import sys from pathlib import Path import numpy as np sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "tier1_scripts")) from _common import import_runner, paired_change, read_endpoints # noqa: E402 CONDITIONS = ("Control", "Hyperosmotic") CASE_DIR_NAMES = {"control": "Control", "hyperosmotic": "Hyperosmotic"} DISPLACEMENT_NAME = "nonlinear_nodal_displacement.npy" CASE_CONFIG_NAME = "case_config.json" # Published natural volume ratio, used for the dry-run preview only; the real # value is read out of the case config the runner builds. DEFAULT_NATURAL_VOLUME_RATIO = 0.7255172413793104 REQUIRED_CODE_CHANGE = """\ Required code change -- `axis_scale_factors` could not be injected ------------------------------------------------------------------ `build_template_case_configs` did not return dicts carrying a "stimulus" block with "volume_ratio", so the runtime wrapper cannot set the natural axis scales. Minimal upstream diff instead of faking the result: run_vaziri_finan_two_group.py, `build_template_case_configs` add a keyword `axis_scale_factors=None` and, when it is not None, write it into the hyperosmotic case's "stimulus" block next to "volume_ratio". run_reynolds_finan_heterogeneous_two_group.py, `run_analysis` add a matching `axis_scale_factors` keyword and forward it to `build_template_case_configs`. Nothing in the solver needs to change: `_target_axis_scales` already returns `cfg.stimulus.axis_scale_factors` when it is set, and `_growth_tensor` already builds `np.diag(target_scales ** exponent)` from it.""" def _lazy_import(name: str, purpose: str): """Import a heavy optional dependency, or explain what is missing and why.""" try: return importlib.import_module(name) except ImportError as exc: raise SystemExit( f"{name} is required for {purpose}, and it is not importable here.\n" "Run this script in the environment that runs the FEM project (the " "one with meshio, scipy and scikit-fem). The solve itself already " "needs all three." ) from exc def _runner(project_root: Path): """Import the upstream two-group runner, or say what the tree is missing.""" try: return import_runner(project_root) except ImportError as exc: raise SystemExit( f"Importing the two-group runner from {project_root} failed: {exc}.\n" "It needs src/nuclear_envelope_fem, run_extreme_finite_deformation.py, " "run_vaziri_finan_two_group.py and validate_literature_reproduction.py, " "none of which ship in the submission zip. Point --project-root at the " "real working tree." ) from exc def axis_scale_factors(volume_ratio: float, anisotropy: float) -> np.ndarray: """Per-axis natural stretches whose product is exactly `volume_ratio`.""" exponents = np.asarray( [(1.0 - anisotropy) / 3.0, (1.0 - anisotropy) / 3.0, (1.0 + 2.0 * anisotropy) / 3.0] ) return np.power(float(volume_ratio), exponents) def _patch_case_configs(runner, anisotropy: float, applied: list[dict]): """Wrap build_template_case_configs so the natural growth becomes anisotropic. Returns the original callable so the caller can restore it between values. """ original = runner.build_template_case_configs def patched(*args, **kwargs): configs = original(*args, **kwargs) for config in configs: if not isinstance(config, dict) or "stimulus" not in config: raise SystemExit(REQUIRED_CODE_CHANGE) stimulus = config["stimulus"] if "volume_ratio" not in stimulus: raise SystemExit(REQUIRED_CODE_CHANGE) volume_ratio = float(stimulus["volume_ratio"]) scales = axis_scale_factors(volume_ratio, anisotropy) if not np.isclose(float(np.prod(scales)), volume_ratio, rtol=1e-12, atol=0.0): raise SystemExit( "Internal error: axis scales do not reproduce the natural volume " f"ratio ({np.prod(scales)!r} vs {volume_ratio!r})." ) stimulus["axis_scale_factors"] = [float(value) for value in scales] applied.append({"volume_ratio": volume_ratio, "axis_scale_factors": scales.tolist()}) return configs runner.build_template_case_configs = patched return original def _case_directories(run_dir: Path, primary_mesh: str, element_order: str) -> dict[str, Path]: """Find the per-condition solve directories inside one run folder.""" wanted = f"{primary_mesh}/{element_order}" found: dict[str, Path] = {} for path in sorted(run_dir.rglob(DISPLACEMENT_NAME)): case_dir = path.parent condition = CASE_DIR_NAMES.get(case_dir.name) if condition is None or wanted not in case_dir.as_posix(): continue found.setdefault(condition, case_dir) return found def _read_imposed_scales(case_dir: Path) -> list[float] | None: config_path = case_dir / CASE_CONFIG_NAME if not config_path.exists(): return None stimulus = json.loads(config_path.read_text(encoding="utf-8")).get("stimulus", {}) scales = stimulus.get("axis_scale_factors") return [float(value) for value in scales] if scales is not None else None def _semi_axes(points: np.ndarray, cells: np.ndarray) -> dict: """Half-extents and inertia-equivalent semi-axes of a tetrahedral body.""" half_extent = 0.5 * (points.max(axis=0) - points.min(axis=0)) tetra = points[cells] volumes = ( np.einsum( "ij,ij->i", tetra[:, 1] - tetra[:, 0], np.cross(tetra[:, 2] - tetra[:, 0], tetra[:, 3] - tetra[:, 0]), ) / 6.0 ) volumes = np.abs(volumes) total = float(volumes.sum()) centroids = tetra.mean(axis=1) centroid = np.einsum("i,ij->j", volumes, centroids) / total # Exact tetrahedron second moment: int x_i x_j dV = V/20 (sum_k v_ki v_kj + 16 c_i c_j). vertex_moment = np.einsum("nki,nkj->nij", tetra, tetra) centroid_moment = 16.0 * np.einsum("ni,nj->nij", centroids, centroids) second_moment = np.einsum("n,nij->ij", volumes / 20.0, vertex_moment + centroid_moment) covariance = second_moment / total - np.outer(centroid, centroid) eigenvalues, eigenvectors = np.linalg.eigh(covariance) equivalent = np.sqrt(np.maximum(5.0 * eigenvalues, 0.0)) # Report principal axes in x/y/z order by dominant alignment, not eigenvalue order. order = np.argmax(np.abs(eigenvectors), axis=0) axes = np.zeros(3) for column, axis in enumerate(order): axes[axis] = equivalent[column] if np.any(axes == 0.0): # Two principal directions claimed the same coordinate axis; fall back to the # nominal ordering, which is valid only because 7.5 > 5.5 > 3.5 um by construction. axes = np.sort(equivalent)[::-1] return { "half_extent_um": half_extent.tolist(), "equivalent_semi_axes_um": axes.tolist(), "volume_um3": total, "centroid_um": centroid.tolist(), } def _measure_axis_ratios(run_dir, mesh_path, primary_mesh, element_order) -> dict: """Deformed semi-axes per condition, from the stored vertex displacements.""" meshio = _lazy_import("meshio", "reading the mesh behind the stored displacements") case_dirs = _case_directories(Path(run_dir), primary_mesh, element_order) if set(case_dirs) != set(CONDITIONS): missing = sorted(set(CONDITIONS) - set(case_dirs)) return {"error": f"no solved case directory for {missing} under {run_dir}"} mesh = meshio.read(mesh_path) points = np.asarray(mesh.points, dtype=float)[:, :3] cells = None for block in mesh.cells: if block.type == "tetra": cells = np.asarray(block.data, dtype=np.int64) if cells is None: return {"error": f"no tetrahedral block in {mesh_path}"} measured: dict = {"mesh": str(mesh_path), "reference": _semi_axes(points, cells)} for condition, case_dir in case_dirs.items(): displacement = np.load(case_dir / DISPLACEMENT_NAME) if displacement.shape != points.shape: return { "error": ( f"{case_dir/DISPLACEMENT_NAME} has shape {displacement.shape}, mesh has " f"{points.shape}; the stored displacement does not belong to this mesh." ) } measured[condition] = { **_semi_axes(points + displacement, cells), "case_dir": str(case_dir), "imposed_axis_scale_factors": _read_imposed_scales(case_dir), } control = np.asarray(measured["Control"]["equivalent_semi_axes_um"]) hyper = np.asarray(measured["Hyperosmotic"]["equivalent_semi_axes_um"]) stretch = hyper / control measured["realized_axis_stretch_hyper_over_control"] = stretch.tolist() measured["control_axis_ratios_zx_zy_yx"] = [ float(control[2] / control[0]), float(control[2] / control[1]), float(control[1] / control[0]), ] measured["hyperosmotic_axis_ratios_zx_zy_yx"] = [ float(hyper[2] / hyper[0]), float(hyper[2] / hyper[1]), float(hyper[1] / hyper[0]), ] measured["realized_height_to_width_ratio_change"] = float(stretch[2] / stretch[0]) return measured def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument( "--project-root", type=Path, required=True, help="Working tree containing src/nuclear_envelope_fem, meshes and results.", ) parser.add_argument("--out-root", type=Path, required=True) parser.add_argument( "--anisotropy", type=float, nargs="+", default=(0.0, 0.25, 0.5, 0.75, 1.0), help="0 reproduces the published isotropic case; 1 puts all shrinkage in z.", ) parser.add_argument( "--allow-extreme", action="store_true", help="Permit a > 1, which imposes in-plane natural expansion.", ) parser.add_argument("--primary-mesh", default="intermediate_h0p75") parser.add_argument("--element-order", choices=("p1", "p2"), default="p1") parser.add_argument("--mesh-path", type=Path, default=None) parser.add_argument( "--natural-volume-ratio", type=float, default=DEFAULT_NATURAL_VOLUME_RATIO, help="Preview value only; the run uses whatever the case config carries.", ) parser.add_argument("--sparse-backend", choices=("cpu", "cuda", "auto"), default="cpu") parser.add_argument( "--rigid-constraint-linear-solver", choices=("augmented_direct", "projected_minres"), default="augmented_direct", ) parser.add_argument("--resume", action="store_true") parser.add_argument("--json-out", type=Path, default=None) parser.add_argument("--markdown-out", type=Path, default=None) parser.add_argument( "--dry-run", action="store_true", help="Print the imposed growth for every value and exit without solving.", ) args = parser.parse_args() extreme = [value for value in args.anisotropy if value > 1.0 or value < 0.0] if extreme and not args.allow_extreme: raise SystemExit( f"anisotropy values {extreme} lie outside [0, 1]. a > 1 imposes in-plane " "natural EXPANSION and a < 0 flattens the wrong way; pass --allow-extreme " "if that is really the intent." ) lines: list[str] = [] def emit(text: str = "") -> None: print(text) lines.append(text) discretization = f"{args.primary_mesh}_{args.element_order.upper()}" preview = float(args.natural_volume_ratio) report: dict = { "project_root": str(args.project_root), "out_root": str(args.out_root), "discretization": discretization, "preview_natural_volume_ratio": preview, "parameterisation": "s = V ** [(1-a)/3, (1-a)/3, (1+2a)/3]; prod(s) = V; s_z/s_x = V**a", "runs": {}, } emit("# Anisotropic natural shrinkage and the measurable axis ratios") emit() emit(f"Discretisation `{discretization}`, natural volume ratio {preview:.6f} " "(preview; the run reads its own).") emit() emit("| a | s_x | s_y | s_z | product (must equal V) | imposed s_z/s_x = V^a | out dir |") emit("|---:|---:|---:|---:|---:|---:|---|") planned = [] for value in args.anisotropy: scales = axis_scale_factors(preview, value) out_dir = args.out_root / f"anisotropy_{value:.2f}".replace(".", "p") planned.append({"anisotropy": value, "out_dir": out_dir, "preview_scales": scales}) emit( f"| {value:.2f} | {scales[0]:.6f} | {scales[1]:.6f} | {scales[2]:.6f} | " f"{float(np.prod(scales)):.6f} | {preview**value:.6f} | `{out_dir.name}` |" ) emit() emit("Control keeps V = 1, so every scale is exactly 1 there and only the hyperosmotic " "natural state moves. `natural_volume_ratio` in the summary CSV is the product of " "these three numbers and must stay at its published value in every row.") if args.dry_run: emit() emit(f"--dry-run: nothing solved. {len(planned)} paired runs would go to " f"`{args.out_root}`, each a Control/Hyperosmotic pair at {discretization}.") if args.json_out: report["planned"] = [ { "anisotropy": item["anisotropy"], "out_dir": str(item["out_dir"]), "preview_axis_scale_factors": item["preview_scales"].tolist(), } for item in planned ] args.json_out.write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8") print(f"\nWrote {args.json_out}", file=sys.stderr) if args.markdown_out: args.markdown_out.write_text("\n".join(lines) + "\n", encoding="utf-8") print(f"Wrote {args.markdown_out}", file=sys.stderr) return 0 runner = _runner(args.project_root) mesh_path = args.mesh_path if mesh_path is None: table = getattr(runner, "MESHES", {}) mesh_path = table.get(args.primary_mesh) if mesh_path is None or not Path(mesh_path).exists(): raise SystemExit( f"Could not locate the mesh for {args.primary_mesh}; pass --mesh-path. " "The axis ratios are reconstructed from the mesh plus the stored vertex " "displacements, so the mesh is required." ) args.out_root.mkdir(parents=True, exist_ok=True) for item in planned: value = item["anisotropy"] out_dir = item["out_dir"] print(f"\n=== anisotropy {value:.2f} -> {out_dir} ===", flush=True) applied: list[dict] = [] original = _patch_case_configs(runner, value, applied) try: runner.run_analysis( out_dir, primary_mesh=args.primary_mesh, run_coarse_p2=False, run_intermediate_p2=(args.element_order == "p2"), resume=args.resume, sparse_backend=args.sparse_backend, rigid_constraint_linear_solver=args.rigid_constraint_linear_solver, ) finally: runner.build_template_case_configs = original entry: dict = {"out_dir": str(out_dir), "applied_axis_scale_factors": applied} try: rows = read_endpoints(out_dir) control, hyper, difference, percent = paired_change(rows, discretization) entry.update( { "control_mn_per_m": control, "hyperosmotic_mn_per_m": hyper, "absolute_difference_mn_per_m": difference, "percent_change": percent, "natural_volume_ratio": float( rows[(discretization, "Hyperosmotic")]["natural_volume_ratio"] ), "realized_outer_volume_ratio": float( rows[(discretization, "Hyperosmotic")]["realized_outer_volume_ratio"] ), "compression_area_fraction": float( rows[(discretization, "Hyperosmotic")]["compression_area_fraction"] ), } ) except (FileNotFoundError, KeyError) as exc: entry["endpoint_error"] = str(exc) entry["geometry"] = _measure_axis_ratios( out_dir, mesh_path, args.primary_mesh, args.element_order ) report["runs"][f"{value:.2f}"] = entry emit() emit("## Mechanical endpoint") emit() emit("| a | Control (mN/m) | Hyperosmotic (mN/m) | difference (mN/m) | change | " "natural V | realised V | compression area |") emit("|---:|---:|---:|---:|---:|---:|---:|---:|") endpoint_errors = [] for key, entry in report["runs"].items(): if "endpoint_error" in entry: emit(f"| {key} | -- | -- | -- | -- | -- | -- | -- |") endpoint_errors.append(f"a = {key}: {entry['endpoint_error']}") continue emit( f"| {key} | {entry['control_mn_per_m']:.6f} | {entry['hyperosmotic_mn_per_m']:.6f} | " f"{entry['absolute_difference_mn_per_m']:+.6f} | {entry['percent_change']:+.2f}% | " f"{entry['natural_volume_ratio']:.6f} | {entry['realized_outer_volume_ratio']:.6f} | " f"{entry['compression_area_fraction']:.3f} |" ) for message in endpoint_errors: emit(f"\n- endpoint unavailable, {message}") volumes = { round(entry["natural_volume_ratio"], 9) for entry in report["runs"].values() if "natural_volume_ratio" in entry } emit() emit( f"- prescribed natural volume ratio across the sweep: {sorted(volumes)} -- one value " "means the sweep changed shape only." if len(volumes) == 1 else f"- WARNING: the natural volume ratio moved across the sweep ({sorted(volumes)}). " "The rows are not comparable; the anisotropy injection changed more than the shape." ) emit() emit("## Measurable prediction: deformed semi-axes and their ratios") emit() emit("| a | Control z/x | Hyper z/x | realised z/x change | imposed V^a | " "realised axis stretch (x, y, z) |") emit("|---:|---:|---:|---:|---:|---|") for key, entry in report["runs"].items(): geometry = entry.get("geometry", {}) if "error" in geometry: emit(f"| {key} | -- | -- | -- | -- | {geometry['error']} |") continue control_ratio = geometry["control_axis_ratios_zx_zy_yx"][0] hyper_ratio = geometry["hyperosmotic_axis_ratios_zx_zy_yx"][0] stretch = geometry["realized_axis_stretch_hyper_over_control"] imposed = ( entry.get("natural_volume_ratio", float("nan")) ** float(key) if "natural_volume_ratio" in entry else float("nan") ) emit( f"| {key} | {control_ratio:.4f} | {hyper_ratio:.4f} | " f"{geometry['realized_height_to_width_ratio_change']:.4f} | {imposed:.4f} | " f"{stretch[0]:.4f}, {stretch[1]:.4f}, {stretch[2]:.4f} |" ) emit() emit("Semi-axes are the inertia-equivalent ellipsoid of the deformed body " "(a_i = sqrt(5 * eigenvalue) of the volume covariance); half-extents are in the " "JSON alongside them. The realised z/x change is what a confocal z-stack measures " "before and after the 580 mOsm step, and it is NOT the imposed V^a: the difference " "between those two columns is the mechanical response the experiment would test.") emit() emit("Report the a = 0 row as the published model, and quote the range of the realised " "z/x change over the sweep as the prediction's sensitivity to an assumption the " "manuscript currently makes silently.") if args.json_out: args.json_out.write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8") print(f"\nWrote {args.json_out}", file=sys.stderr) if args.markdown_out: args.markdown_out.write_text("\n".join(lines) + "\n", encoding="utf-8") print(f"Wrote {args.markdown_out}", file=sys.stderr) return 0 if __name__ == "__main__": raise SystemExit(main())