Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 74 additions & 17 deletions domains/BioMesh/BioMesh_Tests.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<VascularEdge>(resistance, 1.0);
auto edge2 = std::make_shared<VascularEdge>(resistance, 1.0);
auto edge1 = std::make_shared<VascularEdge>(0, 1, capacity, 1.0);
auto edge2 = std::make_shared<VascularEdge>(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");
}
}

Expand Down Expand Up @@ -50,48 +50,105 @@ 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);
CHECK(matrix.payload_eluted, "Payload MUST elute when both conditions are met spatially and temporally.");
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<VascularEdge>(0, 1, 1.0, 1.0);
auto e12 = std::make_shared<VascularEdge>(1, 2, 1.0, 1.0);
auto e23 = std::make_shared<VascularEdge>(2, 3, 1.0, 1.0);
auto e30 = std::make_shared<VascularEdge>(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;
}
64 changes: 63 additions & 1 deletion domains/BioMesh/OsteoMesh.cpp
Original file line number Diff line number Diff line change
@@ -1,15 +1,77 @@
#include "OsteoMesh.hpp"
#include <cmath>
#include <stdexcept>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <algorithm>
#include <random>

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<VascularEdge> edge) {
adjacent_edges.push_back(edge);
}

void OsteoMeshNetwork::add_node(SkeletalNode node) {
nodes.push_back(node);
}

void OsteoMeshNetwork::add_edge(std::shared_ptr<VascularEdge> 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<Eigen::MatrixXd> 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) {
Expand Down
27 changes: 21 additions & 6 deletions domains/BioMesh/OsteoMesh.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::shared_ptr<VascularEdge>> 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<VascularEdge> edge);
};

class OsteoMeshNetwork {
public:
std::vector<SkeletalNode> nodes;
std::vector<std::shared_ptr<VascularEdge>> edges;

void add_node(SkeletalNode node);
void add_edge(std::shared_ptr<VascularEdge> 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
Expand Down
Loading