diff --git a/domains/BioMesh/BioMesh_Tests.cpp b/domains/BioMesh/BioMesh_Tests.cpp index 366a6db..08004a3 100644 --- a/domains/BioMesh/BioMesh_Tests.cpp +++ b/domains/BioMesh/BioMesh_Tests.cpp @@ -11,17 +11,17 @@ } \ } while (false) -void setup_vaso_occlusive_crisis(SkeletalNode& node, double resistance) { +void setup_vaso_occlusive_crisis(SkeletalNode& node, double capacity) { node.adjacent_edges.clear(); - auto edge1 = std::make_shared(resistance, 1.0); - auto edge2 = std::make_shared(resistance, 1.0); + auto edge1 = std::make_shared(0, 1, capacity, 1.0); + auto edge2 = std::make_shared(1, 0, capacity, 1.0); node.add_edge(edge1); node.add_edge(edge2); } -void verify_edges(const SkeletalNode& node, double expected_resistance) { +void verify_edges(const SkeletalNode& node, double expected_capacity) { for (const auto& edge : node.adjacent_edges) { - CHECK(std::abs(edge->flow_resistance - expected_resistance) < 1e-9, "Edge flow resistance does not match expected"); + CHECK(std::abs(edge->flow_capacity - expected_capacity) < 1e-9, "Edge flow capacity does not match expected"); } } @@ -50,41 +50,41 @@ int main() { // Case 1: Neither condition met (normal pH, low protease) { - SkeletalNode node(7.4, 0.5); // pH=7.4, protease=0.5 - setup_vaso_occlusive_crisis(node, 10.0); + SkeletalNode node(0, 7.4, 0.5); // id=0, pH=7.4, protease=0.5 + setup_vaso_occlusive_crisis(node, 0.1); // capacity = 0.1 (degraded) MetaboJointMatrix matrix(beta, w0, E_base, F_enz); matrix.process_tick(node); CHECK(!matrix.payload_eluted, "Payload should NOT elute when neither condition is met."); - verify_edges(node, 10.0); + verify_edges(node, 0.1); } // Case 2: Only Condition 1 met (low pH, low protease) { - SkeletalNode node(6.8, 0.5); // pH=6.8 (acidosis), protease=0.5 - setup_vaso_occlusive_crisis(node, 10.0); + SkeletalNode node(0, 6.8, 0.5); // id=0, pH=6.8 (acidosis), protease=0.5 + setup_vaso_occlusive_crisis(node, 0.1); MetaboJointMatrix matrix(beta, w0, E_base, F_enz); matrix.process_tick(node); CHECK(!matrix.payload_eluted, "Payload should NOT elute when only pH condition is met."); - verify_edges(node, 10.0); + verify_edges(node, 0.1); } // Case 3: Only Condition 2 met (normal pH, high protease) { - SkeletalNode node(7.4, 2.0); // pH=7.4, protease=2.0 (high) - setup_vaso_occlusive_crisis(node, 10.0); + SkeletalNode node(0, 7.4, 2.0); // id=0, pH=7.4, protease=2.0 (high) + setup_vaso_occlusive_crisis(node, 0.1); MetaboJointMatrix matrix(beta, w0, E_base, F_enz); matrix.process_tick(node); CHECK(!matrix.payload_eluted, "Payload should NOT elute when only enzymatic condition is met."); - verify_edges(node, 10.0); + verify_edges(node, 0.1); } // Case 4: Both conditions met (low pH, high protease) { - SkeletalNode node(6.8, 2.0); // pH=6.8 (acidosis), protease=2.0 (high) - setup_vaso_occlusive_crisis(node, 10.0); + SkeletalNode node(0, 6.8, 2.0); // id=0, pH=6.8 (acidosis), protease=2.0 (high) + setup_vaso_occlusive_crisis(node, 0.1); MetaboJointMatrix matrix(beta, w0, E_base, F_enz); matrix.process_tick(node); @@ -92,6 +92,63 @@ int main() { verify_edges(node, 1.0); // baseline is 1.0 } - std::cout << "All Vaso-Occlusive Crisis tests passed successfully. AND gate logic verified." << std::endl; + std::cout << "All AND gate logic tests passed successfully." << std::endl; + + std::cout << "Running Graph-Theoretic Vaso-Occlusive Loop Tests..." << std::endl; + + // Fiedler Value & Vaso-Occlusive Network Degradation / Recovery Test + { + OsteoMeshNetwork network; + + // Create a small ring topology (4 nodes) + SkeletalNode n0(0, 7.4, 0.5); + SkeletalNode n1(1, 7.4, 0.5); + SkeletalNode n2(2, 6.8, 2.0); // n2 is our target niche node (acidosis + high protease) + SkeletalNode n3(3, 7.4, 0.5); + + // Baseline capacity = 1.0 + auto e01 = std::make_shared(0, 1, 1.0, 1.0); + auto e12 = std::make_shared(1, 2, 1.0, 1.0); + auto e23 = std::make_shared(2, 3, 1.0, 1.0); + auto e30 = std::make_shared(3, 0, 1.0, 1.0); + + n0.add_edge(e01); n0.add_edge(e30); + n1.add_edge(e01); n1.add_edge(e12); + n2.add_edge(e12); n2.add_edge(e23); + n3.add_edge(e23); n3.add_edge(e30); + + network.add_node(n0); + network.add_node(n1); + network.add_node(n2); + network.add_node(n3); + + network.add_edge(e01); + network.add_edge(e12); + network.add_edge(e23); + network.add_edge(e30); + + double baseline_fiedler = network.compute_fiedler_value(); + std::cout << "Baseline Fiedler value (λ₂): " << baseline_fiedler << std::endl; + CHECK(baseline_fiedler > 0.0, "Baseline graph should be connected with a positive Fiedler value"); + + // Simulate VOC - randomly degrade edge capacities + network.simulate_vaso_occlusive_loop(20, 0.2); // Degrade 20 times by 0.2 + double degraded_fiedler = network.compute_fiedler_value(); + std::cout << "Degraded Fiedler value (λ₂) after Vaso-Occlusive Crisis: " << degraded_fiedler << std::endl; + CHECK(degraded_fiedler < baseline_fiedler, "Fiedler value should drop after simulated vaso-occlusion"); + + // Trigger the payload on our target node n2 + MetaboJointMatrix payload(beta, w0, E_base, F_enz); + payload.process_tick(n2); // Node 2 meets both conditions (pH 6.8, Protease 2.0) + CHECK(payload.payload_eluted, "Payload should elute at target node n2"); + + // Compute Fiedler value after partial recovery + double recovered_fiedler = network.compute_fiedler_value(); + std::cout << "Recovered Fiedler value (λ₂) after Payload Elution: " << recovered_fiedler << std::endl; + CHECK(recovered_fiedler > degraded_fiedler, "Fiedler value should rise after restorative payload elution"); + } + + std::cout << "Graph Laplacian and Network Recovery tests passed successfully." << std::endl; + return 0; } diff --git a/domains/BioMesh/OsteoMesh.cpp b/domains/BioMesh/OsteoMesh.cpp index 6979d3d..e6bb06c 100644 --- a/domains/BioMesh/OsteoMesh.cpp +++ b/domains/BioMesh/OsteoMesh.cpp @@ -1,15 +1,77 @@ #include "OsteoMesh.hpp" #include #include +#include +#include +#include +#include void VascularEdge::reset_to_baseline() { - flow_resistance = baseline_resistance; + flow_capacity = baseline_capacity; +} + +void VascularEdge::degrade(double amount) { + flow_capacity -= amount; + if (flow_capacity < 0.0) { + flow_capacity = 0.0; + } } void SkeletalNode::add_edge(std::shared_ptr edge) { adjacent_edges.push_back(edge); } +void OsteoMeshNetwork::add_node(SkeletalNode node) { + nodes.push_back(node); +} + +void OsteoMeshNetwork::add_edge(std::shared_ptr edge) { + edges.push_back(edge); +} + +double OsteoMeshNetwork::compute_fiedler_value() { + int n = nodes.size(); + if (n == 0) return 0.0; + + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n, n); + Eigen::MatrixXd D = Eigen::MatrixXd::Zero(n, n); + + for (const auto& edge : edges) { + int u = edge->u; + int v = edge->v; + A(u, v) += edge->flow_capacity; + A(v, u) += edge->flow_capacity; + } + + for (int i = 0; i < n; ++i) { + double degree = 0.0; + for (int j = 0; j < n; ++j) { + degree += A(i, j); + } + D(i, i) = degree; + } + + Eigen::MatrixXd L = D - A; + + Eigen::SelfAdjointEigenSolver solver(L); + Eigen::VectorXd eigenvalues = solver.eigenvalues(); + + if (n < 2) return 0.0; + return eigenvalues(1); // Fiedler value is the second smallest eigenvalue +} + +void OsteoMeshNetwork::simulate_vaso_occlusive_loop(int iterations, double degrade_amount) { + std::random_device rd; + std::mt19937 gen(rd()); + if (edges.empty()) return; + std::uniform_int_distribution<> distrib(0, edges.size() - 1); + + for (int i = 0; i < iterations; ++i) { + int rand_edge_idx = distrib(gen); + edges[rand_edge_idx]->degrade(degrade_amount); + } +} + MetaboJointMatrix::MetaboJointMatrix(double b, double w0, double E_base, double F_enz) : beta(b), omega0(w0), E_baseline(E_base), F_enzyme_impulse(F_enz), pH_threshold(7.2), protease_threshold(1.0), payload_eluted(false) { diff --git a/domains/BioMesh/OsteoMesh.hpp b/domains/BioMesh/OsteoMesh.hpp index a2c368d..5495436 100644 --- a/domains/BioMesh/OsteoMesh.hpp +++ b/domains/BioMesh/OsteoMesh.hpp @@ -5,27 +5,42 @@ class VascularEdge { public: - double flow_resistance; - double baseline_resistance; + int u; + int v; + double flow_capacity; + double baseline_capacity; - VascularEdge(double resistance = 1.0, double baseline = 1.0) - : flow_resistance(resistance), baseline_resistance(baseline) {} + VascularEdge(int u_id = 0, int v_id = 0, double capacity = 1.0, double baseline = 1.0) + : u(u_id), v(v_id), flow_capacity(capacity), baseline_capacity(baseline) {} void reset_to_baseline(); + void degrade(double amount); }; class SkeletalNode { public: + int id; double local_pH; double local_protease_concentration; std::vector> adjacent_edges; - SkeletalNode(double pH = 7.4, double protease = 0.0) - : local_pH(pH), local_protease_concentration(protease) {} + SkeletalNode(int node_id = 0, double pH = 7.4, double protease = 0.0) + : id(node_id), local_pH(pH), local_protease_concentration(protease) {} void add_edge(std::shared_ptr edge); }; +class OsteoMeshNetwork { +public: + std::vector nodes; + std::vector> edges; + + void add_node(SkeletalNode node); + void add_edge(std::shared_ptr edge); + double compute_fiedler_value(); + void simulate_vaso_occlusive_loop(int iterations, double degrade_amount); +}; + class MetaboJointMatrix { public: double beta; // Must be < 0 for bistable regime