Conserve biochemical biomass across growth, division, and mechanics

cpp/core/coupled_rates.cpp

@@ -96,6 +96,7 @@ void CoupledRatePlan::validate() const {
       case RateOp::growth_rate:
       case RateOp::cell_type:
       case RateOp::cell_volume:
+      case RateOp::cell_volume_change_rate:
       case RateOp::cell_surface_area:
         break;
       case RateOp::negate:

cpp/core/mechanics_integration.cpp

@@ -72,12 +72,14 @@ void integrate_mechanics_result(WorldState& state, const MechanicsSolveResult& r
       throw std::invalid_argument("desired length increments must be finite and non-negative");
     }
     const auto applied_length_increment =
-        cell.fixed ? desired_increment : std::max(0.0F, desired_increment + correction.length);
+        desired_length_increments.empty()
+            ? 0.0F
+            : (cell.fixed ? desired_increment
+                          : std::max(0.0F, desired_increment + correction.length));
     const auto new_position = cell.fixed ? cell.position : cell.position + correction.translation;
-    const auto new_direction =
-        cell.fixed ? cell.direction
-                   : rotate_axis_angle(cell.direction, correction.rotation,
-                                       parameters.max_rotation_radians);
+    const auto new_direction = cell.fixed ? cell.direction
+                                          : rotate_axis_angle(cell.direction, correction.rotation,
+                                                              parameters.max_rotation_radians);
     const auto new_length = cell.length + applied_length_increment;
     if (!finite(new_position) || !finite(new_direction) || !std::isfinite(new_length)) {
       throw std::overflow_error("mechanics integration produced non-finite geometry");

cpp/core/species.cpp

@@ -85,6 +85,7 @@ void SpeciesRatePlan::validate() const {
       case RateOp::growth_rate:
       case RateOp::cell_type:
       case RateOp::cell_volume:
+      case RateOp::cell_volume_change_rate:
       case RateOp::cell_surface_area:
         break;
       case RateOp::negate:
@@ -125,6 +126,8 @@ void SpeciesRatePlan::validate() const {
 }
 
 float effective_cell_volume(float length, float radius) noexcept {
+  // Conserved biomass volume for the endpoint-preserving division rule.
+  // This is a biochemical measure, not the geometric capsule volume.
   constexpr float pi = 3.14159265358979323846F;
   return pi * radius * radius * (length + 2.0F * radius);
 }

cpp/cpu/cpu_coupled.cpp

@@ -16,7 +16,7 @@ namespace {
 float evaluate_instruction(const RateInstruction& instruction, std::span<const float> workspace,
                            std::span<const float> species, std::span<const float> signals,
                            const CellGeometryView& geometry, const CellAttributeView& attributes,
-                           std::size_t cell) {
+                           std::size_t cell, float volume_change_rate) {
   switch (instruction.operation) {
     case RateOp::constant:
       return instruction.value;
@@ -38,6 +38,8 @@ float evaluate_instruction(const RateInstruction& instruction, std::span<const f
       return attributes.growth_rates[cell];
     case RateOp::cell_type:
       return static_cast<float>(attributes.cell_types[cell]);
+    case RateOp::cell_volume_change_rate:
+      return volume_change_rate;
     case RateOp::cell_volume:
       return effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]);
     case RateOp::cell_surface_area:
@@ -161,8 +163,13 @@ SignalSolveReport advance_coupled_cpu(WorldState& state, SignalGrid& grid,
     const auto cell_signals =
         std::span<const float>(sampled).subspan(cell * plan.signal_count(), plan.signal_count());
     for (std::size_t index = 0; index < plan.instructions().size(); ++index) {
-      workspace[index] = evaluate_instruction(plan.instructions()[index], workspace, cell_species,
-                                              cell_signals, geometry, attributes, cell);
+      workspace[index] = evaluate_instruction(
+          plan.instructions()[index], workspace, cell_species, cell_signals, geometry, attributes,
+          cell,
+          dt == 0.0F ? 0.0F
+                     : (effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]) -
+                        effective_cell_volume(previous_lengths[cell], geometry.radii[cell])) /
+                           dt);
       if (!std::isfinite(workspace[index])) {
         throw std::domain_error("coupled rate instruction " + std::to_string(index) +
                                 " produced a non-finite value");

cpp/cpu/cpu_species.cpp

@@ -13,7 +13,8 @@ namespace {
 
 float evaluate_instruction(const RateInstruction& instruction, std::span<const float> workspace,
                            std::span<const float> species, const CellGeometryView& geometry,
-                           const CellAttributeView& attributes, std::size_t cell) {
+                           const CellAttributeView& attributes, std::size_t cell,
+                           float volume_change_rate) {
   switch (instruction.operation) {
     case RateOp::constant:
       return instruction.value;
@@ -35,6 +36,8 @@ float evaluate_instruction(const RateInstruction& instruction, std::span<const f
       return attributes.growth_rates[cell];
     case RateOp::cell_type:
       return static_cast<float>(attributes.cell_types[cell]);
+    case RateOp::cell_volume_change_rate:
+      return volume_change_rate;
     case RateOp::cell_volume:
       return effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]);
     case RateOp::cell_surface_area:
@@ -119,8 +122,12 @@ void advance_species_cpu(WorldState& state, const SpeciesRatePlan& plan,
     const auto cell_species =
         std::span<const float>(next_levels).subspan(offset, state.species_count());
     for (std::size_t index = 0; index < plan.instructions().size(); ++index) {
-      workspace[index] = evaluate_instruction(plan.instructions()[index], workspace, cell_species,
-                                              geometry, attributes, cell);
+      workspace[index] = evaluate_instruction(
+          plan.instructions()[index], workspace, cell_species, geometry, attributes, cell,
+          dt == 0.0F ? 0.0F
+                     : (effective_cell_volume(geometry.lengths[cell], geometry.radii[cell]) -
+                        effective_cell_volume(previous_lengths[cell], geometry.radii[cell])) /
+                           dt);
       if (!std::isfinite(workspace[index])) {
         throw std::domain_error("species rate instruction " + std::to_string(index) +
                                 " produced a non-finite value");

cpp/cuda/kernels/coupled_rates.cu

@@ -92,7 +92,8 @@ __device__ float sample_signal(const float* levels, SignalGridShapeGpu shape, fl
 
 __device__ float evaluate_instruction(const RateInstructionGpu& instruction, const float* workspace,
                                       const float* species, const float* signals, float4 center,
-                                      float4 geometry, float growth_rate, std::int32_t cell_type) {
+                                      float4 geometry, float growth_rate, std::int32_t cell_type,
+                                      float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -112,6 +113,8 @@ __device__ float evaluate_instruction(const RateInstructionGpu& instruction, con
       return growth_rate;
     case 8:
       return static_cast<float>(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -190,7 +193,11 @@ __global__ void advance_coupled_cells(
   for (std::uint32_t index = 0; index < instruction_count; ++index) {
     const auto value =
         evaluate_instruction(instructions[index], cell_workspace, cell_species, cell_signals,
-                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell]);
+                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell],
+                             dt == 0.0f ? 0.0f
+                                        : (effective_volume(geometry[cell].x, radius) -
+                                           effective_volume(previous_lengths[cell], radius)) /
+                                              dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomicOr(error, 1U);

cpp/cuda/kernels/species.cu

@@ -16,7 +16,8 @@ __device__ float effective_surface_area(float length, float radius) {
 
 __device__ float evaluate_instruction(const RateInstructionGpu& instruction, const float* workspace,
                                       const float* species, float4 center, float4 geometry,
-                                      float growth_rate, std::int32_t cell_type) {
+                                      float growth_rate, std::int32_t cell_type,
+                                      float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -36,6 +37,8 @@ __device__ float evaluate_instruction(const RateInstructionGpu& instruction, con
       return growth_rate;
     case 8:
       return static_cast<float>(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -104,7 +107,11 @@ __global__ void advance_species(float* levels, const float* previous_lengths, co
   for (std::uint32_t index = 0; index < instruction_count; ++index) {
     const auto value =
         evaluate_instruction(instructions[index], cell_workspace, cell_species, centers[cell],
-                             geometry[cell], growth_rates[cell], cell_types[cell]);
+                             geometry[cell], growth_rates[cell], cell_types[cell],
+                             dt == 0.0f ? 0.0f
+                                        : (effective_volume(geometry[cell].x, radius) -
+                                           effective_volume(previous_lengths[cell], radius)) /
+                                              dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomicOr(error, 1U);

cpp/include/cm/species.hpp

@@ -38,6 +38,7 @@ enum class RateOp : std::uint8_t {
   equal = 25,
   select = 26,
   signal = 27,
+  cell_volume_change_rate = 28,
 };
 
 struct RateInstruction {

cpp/metal/kernels/coupled_rates.metal

@@ -99,7 +99,7 @@ float sample_signal(device const float* levels, GridShape shape, float4 origin,
 
 float evaluate_instruction(const RateInstruction instruction, device const float* workspace,
                            device const float* species, device const float* signals, float4 center,
-                           float4 geometry, float growth_rate, int cell_type) {
+                           float4 geometry, float growth_rate, int cell_type, float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -119,6 +119,8 @@ float evaluate_instruction(const RateInstruction instruction, device const float
       return growth_rate;
     case 8:
       return float(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -201,7 +203,9 @@ kernel void advance_coupled_cells(
   for (uint index = 0; index < instruction_count; ++index) {
     float value =
         evaluate_instruction(instructions[index], cell_workspace, cell_species, cell_signals,
-                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell]);
+                             centers[cell], geometry[cell], growth_rates[cell], cell_types[cell],
+            dt == 0.0f ? 0.0f : (effective_volume(geometry[cell].x, radius) -
+                                 effective_volume(previous_lengths[cell], radius)) / dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomic_fetch_or_explicit(error, 1u, memory_order_relaxed);

cpp/metal/kernels/species.metal

@@ -26,7 +26,7 @@ float evaluate_instruction(const RateInstruction instruction,
                            float4 center,
                            float4 geometry,
                            float growth_rate,
-                           int cell_type) {
+                           int cell_type, float volume_change_rate) {
   switch (instruction.operation) {
     case 0:
       return instruction.value;
@@ -46,6 +46,8 @@ float evaluate_instruction(const RateInstruction instruction,
       return growth_rate;
     case 8:
       return float(cell_type);
+    case 28:
+      return volume_change_rate;
     case 9:
       return effective_volume(geometry.x, geometry.y);
     case 10:
@@ -122,7 +124,9 @@ kernel void advance_species(
   for (uint index = 0; index < instruction_count; ++index) {
     float value = evaluate_instruction(instructions[index], cell_workspace, cell_species,
                                        centers[cell], geometry[cell], growth_rates[cell],
-                                       cell_types[cell]);
+                                       cell_types[cell],
+            dt == 0.0f ? 0.0f : (effective_volume(geometry[cell].x, radius) -
+                                 effective_volume(previous_lengths[cell], radius)) / dt);
     cell_workspace[index] = value;
     if (!isfinite(value)) {
       atomic_fetch_or_explicit(error, 1u, memory_order_relaxed);

cpp/python/bindings.cpp

@@ -49,6 +49,7 @@ NB_MODULE(_core, module) {
       .value("GROWTH_RATE", cm::RateOp::growth_rate)
       .value("CELL_TYPE", cm::RateOp::cell_type)
       .value("CELL_VOLUME", cm::RateOp::cell_volume)
+      .value("CELL_VOLUME_CHANGE_RATE", cm::RateOp::cell_volume_change_rate)
       .value("CELL_SURFACE_AREA", cm::RateOp::cell_surface_area)
       .value("ADD", cm::RateOp::add)
       .value("SUBTRACT", cm::RateOp::subtract)
@@ -215,7 +216,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("instructions",
                    [](const cm::SpeciesRatePlan& plan) {
                      return std::vector<cm::RateInstruction>(plan.instructions().begin(),
-                                                              plan.instructions().end());
+                                                             plan.instructions().end());
                    })
       .def_prop_ro("outputs",
                    [](const cm::SpeciesRatePlan& plan) {
@@ -234,7 +235,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("instructions",
                    [](const cm::CoupledRatePlan& plan) {
                      return std::vector<cm::RateInstruction>(plan.instructions().begin(),
-                                                              plan.instructions().end());
+                                                             plan.instructions().end());
                    })
       .def_prop_ro("species_outputs",
                    [](const cm::CoupledRatePlan& plan) {
@@ -271,7 +272,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("contacts",
                    [](const cm::ContactGraph& graph) {
                      return std::vector<cm::CellContact>(graph.contacts().begin(),
-                                                          graph.contacts().end());
+                                                         graph.contacts().end());
                    })
       .def("__len__", &cm::ContactGraph::size)
       .def(
@@ -357,7 +358,7 @@ NB_MODULE(_core, module) {
       .def_prop_ro("contacts",
                    [](const cm::ExternalContactGraph& graph) {
                      return std::vector<cm::ExternalContact>(graph.contacts().begin(),
-                                                              graph.contacts().end());
+                                                             graph.contacts().end());
                    })
       .def("__len__", &cm::ExternalContactGraph::size)
       .def(
@@ -400,8 +401,8 @@ NB_MODULE(_core, module) {
       .def(nb::init<cm::BackendKind, std::size_t, std::size_t, std::uint32_t>(),
            "backend"_a = cm::BackendKind::cpu, "reserved_capacity"_a = 0, "species_count"_a = 0,
            "device_index"_a = 0)
-      .def(nb::init<cm::BackendKind, const cm::SimulationCheckpoint&, std::uint32_t>(),
-           "backend"_a, "checkpoint"_a, "device_index"_a = 0)
+      .def(nb::init<cm::BackendKind, const cm::SimulationCheckpoint&, std::uint32_t>(), "backend"_a,
+           "checkpoint"_a, "device_index"_a = 0)
       .def_prop_ro("backend_info", &cm::Simulation::backend_info)
       .def("supports", &cm::Simulation::supports, "feature"_a)
       .def_prop_ro("time", &cm::Simulation::time)

docs/architecture/0024-biomass-accounting.md

@@ -0,0 +1,11 @@
+# Biomass accounting
+
+The conserved biochemical volume is B = pi r²(l + 2r), with rod centerline length l and radius r. At constant biomass density rho, biomass mass is rho B and an intracellular concentration c represents amount c B. This effective volume is distinct from the geometric capsule volume pi r²l + 4pi r³/3. The latter changes when a parent capsule is replaced by two rounded daughters, so it must not be used interchangeably with B in biomass-dependent feedback.
+
+Native division preserves the parent's outer endpoints and satisfies l1 + l2 = l - 2r. It therefore conserves B exactly in real arithmetic, including unequal division. The division fraction partitions the available cylindrical length, not total biomass. Both daughters inherit concentrations, so intracellular amount is conserved too. Geometry remains a coarse representation of division rather than a volume-preserving model of septum formation.
+
+Growth retains the elongation law dl/dt = g l. Consequently dB/dt = pi r² g l, not g B. `RatePlanBuilder.cell_volume_change_rate()` returns the realized change (B_after - B_before)/dt during the current growth step, with zero at dt=0. A nutrient sink `-rates.cell_volume_change_rate()/yield` therefore consumes exactly the biomass increment divided by yield, up to floating-point error. `rates.cell_volume()` and `microsimulator.biomass.biomass_volume` expose B. A negative concentration is a rejected step, never an invitation to clip away a mass-balance error.
+
+Contact relaxation without explicitly prescribed length increments preserves each rod's length. Numerical overlap corrections therefore cannot manufacture biomass after the biological update. Direct geometry edits and explicitly prescribed mechanical length increments are externally imposed changes; callers must account for their biomass and intracellular amounts.
+
+Colony resistance uses a calibrated biomass-density closure derived from B. Its density must be deposited conservatively over a physical averaging scale independent of voxel size. It is not an exact solid-occupancy reconstruction of individual capsules.

docs/architecture/README.md

@@ -29,3 +29,5 @@ The [numerical contract](numerical-contract.md) is the best starting point for w
 - [Versioned columnar analysis datasets](0013-analysis-datasets.md)
 
 Current backend status and the tests required to support it are documented in [testing and validation](../development/validation.md).
+
+- [Biomass accounting](0024-biomass-accounting.md)

python/src/microsimulator/_core.pyi

@@ -35,6 +35,7 @@ class RateOp(Enum):
     GROWTH_RATE: RateOp
     CELL_TYPE: RateOp
     CELL_VOLUME: RateOp
+    CELL_VOLUME_CHANGE_RATE: RateOp
     CELL_SURFACE_AREA: RateOp
     ADD: RateOp
     SUBTRACT: RateOp

python/src/microsimulator/biomass.py

@@ -0,0 +1,19 @@
+"""The conserved biochemical volume and its distinct geometric counterpart."""
+
+import math
+
+
+def biomass_volume(length: float, radius: float) -> float:
+    """Return pi*r**2*(length + 2*r), conserved by native cell division.
+
+    At constant biomass density, multiplying this measure by density gives
+    biomass mass. It is the same volume used by native concentration dilution.
+    """
+    if not math.isfinite(length) or length < 0 or not math.isfinite(radius) or radius <= 0:
+        raise ValueError("length must be finite and nonnegative; radius finite and positive")
+    return math.pi * radius * radius * (length + 2 * radius)
+
+
+def capsule_volume(length: float, radius: float) -> float:
+    """Return the geometric capsule volume, which is not the biomass measure."""
+    return biomass_volume(length, radius) - (2 / 3) * math.pi * radius**3

python/src/microsimulator/checkpoint.py

@@ -93,6 +93,7 @@ class CheckpointBundle:
     RateOp.GROWTH_RATE: "growth_rate",
     RateOp.CELL_TYPE: "cell_type",
     RateOp.CELL_VOLUME: "cell_volume",
+    RateOp.CELL_VOLUME_CHANGE_RATE: "cell_volume_change_rate",
     RateOp.CELL_SURFACE_AREA: "cell_surface_area",
     RateOp.ADD: "add",
     RateOp.SUBTRACT: "subtract",