// ═══════════════════════════════════════════════════════════════════════════════
// NEUROLOGICAL MATHEMATICS FRAMEWORK v1.0
// Complete Brain-Spine-Neuron Mathematical Foundation
// ═══════════════════════════════════════════════════════════════════════════════
//
// This module implements the full neurological mathematical framework:
//
// LEVEL 1: NEURON (Single Mind = Single Neuron)
//   - Hodgkin-Huxley equations for action potentials
//   - Nernst equation for ion equilibrium
//   - Goldman-Hodgkin-Katz equation for membrane potential
//   - Synaptic transmission (AMPA, NMDA, GABA)
//   - Spike-timing dependent plasticity (STDP)
//
// LEVEL 2: NEURAL NETWORK (Mind Farm = Neural Tissue)
//   - Connectivity patterns (small-world, scale-free)
//   - Population dynamics (Wilson-Cowan equations)
//   - Neural field theory (Amari equations)
//   - Synchronization (Kuramoto model)
//
// LEVEL 3: BRAIN SYSTEM (Full Nervous System)
//   - Cortical columns and layers
//   - Thalamo-cortical loops
//   - Basal ganglia circuits
//   - Cerebellar computation
//   - Spinal cord reflexes
//
// ═══════════════════════════════════════════════════════════════════════════════

#pragma once

#include <cstdint>
#include <cmath>
#include <vector>
#include <array>
#include <functional>

#include "constants.hpp"
#include "fast_math.hpp"

namespace mind_farm {
namespace neuro {

// ═══════════════════════════════════════════════════════════════════════════════
// PHYSICAL CONSTANTS
// ═══════════════════════════════════════════════════════════════════════════════

constexpr float FARADAY = 96485.0f;      // C/mol - Faraday constant
constexpr float GAS_R = 8.314f;          // J/(mol·K) - Gas constant
constexpr float TEMP_K = 310.0f;         // K - Body temperature (37°C)
constexpr float RT_F = GAS_R * TEMP_K / FARADAY;  // ~26.7 mV at 37°C

// Ion concentrations (mM)
constexpr float NA_OUT = 145.0f;         // [Na+]out
constexpr float NA_IN = 12.0f;           // [Na+]in
constexpr float K_OUT = 4.0f;            // [K+]out
constexpr float K_IN = 155.0f;           // [K+]in
constexpr float CL_OUT = 120.0f;         // [Cl-]out
constexpr float CL_IN = 4.0f;            // [Cl-]in
constexpr float CA_OUT = 2.0f;           // [Ca2+]out
constexpr float CA_IN = 0.0001f;         // [Ca2+]in (~100 nM)

// ═══════════════════════════════════════════════════════════════════════════════
// LEVEL 1: SINGLE NEURON MATHEMATICS
// ═══════════════════════════════════════════════════════════════════════════════

// Nernst Equation: E_ion = (RT/zF) * ln([ion]_out/[ion]_in)
// Calculates equilibrium potential for a single ion
inline float nernst(float z, float conc_out, float conc_in) {
    return (RT_F / z) * std::log(conc_out / conc_in) * 1000.0f;  // Convert to mV
}

// Goldman-Hodgkin-Katz Equation
// V_m = (RT/F) * ln((P_Na[Na]_o + P_K[K]_o + P_Cl[Cl]_i) / (P_Na[Na]_i + P_K[K]_i + P_Cl[Cl]_o))
struct GHKEquation {
    float P_Na = 1.0f;   // Na+ permeability (relative)
    float P_K = 20.0f;   // K+ permeability (at rest, ~20x Na+)
    float P_Cl = 5.0f;   // Cl- permeability

    float membrane_potential() const {
        float num = P_Na * NA_OUT + P_K * K_OUT + P_Cl * CL_IN;
        float den = P_Na * NA_IN + P_K * K_IN + P_Cl * CL_OUT;
        return RT_F * std::log(num / den) * 1000.0f;  // mV
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// HODGKIN-HUXLEY MODEL
// The fundamental equations of neuronal action potentials
//
// C_m * dV/dt = I_ext - g_Na*m³*h*(V-E_Na) - g_K*n⁴*(V-E_K) - g_L*(V-E_L)
//
// Gating variables:
//   dx/dt = α_x(V)*(1-x) - β_x(V)*x   for x ∈ {m, h, n}
//
// ═══════════════════════════════════════════════════════════════════════════════

struct HodgkinHuxley {
    // Membrane properties
    float C_m = 1.0f;        // μF/cm² - Membrane capacitance
    float g_Na = 120.0f;     // mS/cm² - Sodium conductance
    float g_K = 36.0f;       // mS/cm² - Potassium conductance
    float g_L = 0.3f;        // mS/cm² - Leak conductance

    // Reversal potentials (mV)
    float E_Na = 50.0f;
    float E_K = -77.0f;
    float E_L = -54.4f;

    // State variables
    float V = -65.0f;        // mV - Membrane potential
    float m = 0.05f;         // Na+ activation
    float h = 0.6f;          // Na+ inactivation
    float n = 0.32f;         // K+ activation

    // Rate functions for gating variables
    // α_m(V) = 0.1(V+40)/(1-exp(-(V+40)/10))
    float alpha_m(float V_) const {
        float x = V_ + 40.0f;
        if (std::abs(x) < 1e-6f) return 1.0f;  // L'Hôpital limit
        return 0.1f * x / (1.0f - std::exp(-x / 10.0f));
    }

    // β_m(V) = 4*exp(-(V+65)/18)
    float beta_m(float V_) const {
        return 4.0f * std::exp(-(V_ + 65.0f) / 18.0f);
    }

    // α_h(V) = 0.07*exp(-(V+65)/20)
    float alpha_h(float V_) const {
        return 0.07f * std::exp(-(V_ + 65.0f) / 20.0f);
    }

    // β_h(V) = 1/(1+exp(-(V+35)/10))
    float beta_h(float V_) const {
        return 1.0f / (1.0f + std::exp(-(V_ + 35.0f) / 10.0f));
    }

    // α_n(V) = 0.01(V+55)/(1-exp(-(V+55)/10))
    float alpha_n(float V_) const {
        float x = V_ + 55.0f;
        if (std::abs(x) < 1e-6f) return 0.1f;
        return 0.01f * x / (1.0f - std::exp(-x / 10.0f));
    }

    // β_n(V) = 0.125*exp(-(V+65)/80)
    float beta_n(float V_) const {
        return 0.125f * std::exp(-(V_ + 65.0f) / 80.0f);
    }

    // Ionic currents
    float I_Na() const { return g_Na * m*m*m * h * (V - E_Na); }
    float I_K() const { return g_K * n*n*n*n * (V - E_K); }
    float I_L() const { return g_L * (V - E_L); }

    // Time step using Euler method
    // dt in ms, I_ext in μA/cm²
    void step(float dt, float I_ext) {
        // Gating variable derivatives
        float dm = alpha_m(V) * (1.0f - m) - beta_m(V) * m;
        float dh = alpha_h(V) * (1.0f - h) - beta_h(V) * h;
        float dn = alpha_n(V) * (1.0f - n) - beta_n(V) * n;

        // Membrane potential derivative
        float dV = (I_ext - I_Na() - I_K() - I_L()) / C_m;

        // Euler integration
        V += dV * dt;
        m += dm * dt;
        h += dh * dt;
        n += dn * dt;

        // Clamp gating variables to [0,1]
        m = std::clamp(m, 0.0f, 1.0f);
        h = std::clamp(h, 0.0f, 1.0f);
        n = std::clamp(n, 0.0f, 1.0f);
    }

    // Check if neuron is spiking (V > 0 mV)
    bool is_spiking() const { return V > 0.0f; }
};

// ═══════════════════════════════════════════════════════════════════════════════
// SYNAPTIC TRANSMISSION
// Models for AMPA, NMDA, and GABA receptors
// ═══════════════════════════════════════════════════════════════════════════════

// AMPA Receptor (fast excitatory)
// I_AMPA = g_AMPA * s * (V - E_AMPA)
// ds/dt = -s/τ_AMPA + Σ_j δ(t - t_j)
struct AMPASynapse {
    float g_max = 0.35f;     // nS - Maximum conductance
    float E_rev = 0.0f;      // mV - Reversal potential
    float tau = 2.0f;        // ms - Decay time constant
    float s = 0.0f;          // Synaptic variable

    float current(float V_post) const {
        return g_max * s * (V_post - E_rev);
    }

    void receive_spike() {
        s += 1.0f;
    }

    void step(float dt) {
        s *= std::exp(-dt / tau);
    }
};

// NMDA Receptor (slow excitatory, voltage-dependent Mg2+ block)
// I_NMDA = g_NMDA * s * B(V) * (V - E_NMDA)
// B(V) = 1 / (1 + [Mg2+]_o * exp(-0.062*V) / 3.57)
struct NMDASynapse {
    float g_max = 0.5f;      // nS
    float E_rev = 0.0f;      // mV
    float tau_rise = 2.0f;   // ms
    float tau_decay = 100.0f;// ms - NMDA is slow!
    float Mg_o = 1.0f;       // mM - External Mg2+ concentration
    float s = 0.0f;
    float x = 0.0f;          // Rise variable

    // Mg2+ block factor
    float B(float V) const {
        return 1.0f / (1.0f + Mg_o * std::exp(-0.062f * V) / 3.57f);
    }

    float current(float V_post) const {
        return g_max * s * B(V_post) * (V_post - E_rev);
    }

    void receive_spike() {
        x += 1.0f;
    }

    void step(float dt) {
        float dx = -x / tau_rise;
        float ds = -s / tau_decay + x / tau_rise;
        x += dx * dt;
        s += ds * dt;
    }
};

// GABA-A Receptor (fast inhibitory)
struct GABAASynapse {
    float g_max = 0.8f;      // nS
    float E_rev = -80.0f;    // mV (Cl- reversal)
    float tau = 5.0f;        // ms
    float s = 0.0f;

    float current(float V_post) const {
        return g_max * s * (V_post - E_rev);
    }

    void receive_spike() {
        s += 1.0f;
    }

    void step(float dt) {
        s *= std::exp(-dt / tau);
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// SPIKE-TIMING DEPENDENT PLASTICITY (STDP)
// Δw = A+ * exp(-Δt/τ+)  if Δt > 0 (post after pre → LTP)
// Δw = -A- * exp(Δt/τ-)  if Δt < 0 (pre after post → LTD)
// ═══════════════════════════════════════════════════════════════════════════════

struct STDP {
    float A_plus = 0.005f;   // LTP amplitude
    float A_minus = 0.00525f;// LTD amplitude (slightly larger for stability)
    float tau_plus = 20.0f;  // ms - LTP time constant
    float tau_minus = 20.0f; // ms - LTD time constant
    float w_max = 1.0f;      // Maximum weight
    float w_min = 0.0f;      // Minimum weight

    // Calculate weight change given spike time difference
    // delta_t = t_post - t_pre
    float delta_w(float delta_t, float w) const {
        if (delta_t > 0) {
            // Post after pre → LTP
            return A_plus * (w_max - w) * std::exp(-delta_t / tau_plus);
        } else {
            // Pre after post → LTD
            return -A_minus * (w - w_min) * std::exp(delta_t / tau_minus);
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// LEVEL 2: NEURAL NETWORK / POPULATION MATHEMATICS
// ═══════════════════════════════════════════════════════════════════════════════

// Wilson-Cowan Equations (population dynamics)
// τ_E * dE/dt = -E + (1 - r_E*E) * S(w_EE*E - w_EI*I + I_E)
// τ_I * dI/dt = -I + (1 - r_I*I) * S(w_IE*E - w_II*I + I_I)
// where S(x) = 1/(1 + exp(-a*(x-θ))) is the sigmoid
struct WilsonCowan {
    // Time constants (ms)
    float tau_E = 10.0f;
    float tau_I = 20.0f;

    // Connection weights
    float w_EE = 12.0f;   // E → E
    float w_EI = 4.0f;    // I → E (inhibition of excitatory)
    float w_IE = 13.0f;   // E → I
    float w_II = 11.0f;   // I → I

    // Refractory parameters
    float r_E = 1.0f;
    float r_I = 1.0f;

    // Sigmoid parameters
    float a = 1.0f;
    float theta = 4.0f;

    // External inputs
    float I_E = 0.0f;
    float I_I = 0.0f;

    // Population activities
    float E = 0.1f;       // Excitatory population
    float I = 0.1f;       // Inhibitory population

    float sigmoid(float x) const {
        return 1.0f / (1.0f + std::exp(-a * (x - theta)));
    }

    void step(float dt) {
        float input_E = w_EE * E - w_EI * I + I_E;
        float input_I = w_IE * E - w_II * I + I_I;

        float dE = (-E + (1.0f - r_E * E) * sigmoid(input_E)) / tau_E;
        float dI = (-I + (1.0f - r_I * I) * sigmoid(input_I)) / tau_I;

        E += dE * dt;
        I += dI * dt;

        E = std::clamp(E, 0.0f, 1.0f);
        I = std::clamp(I, 0.0f, 1.0f);
    }
};

// Neural Field Theory (Amari Equation)
// τ * ∂u(x,t)/∂t = -u(x,t) + ∫w(x-x')*f(u(x',t))dx' + s(x,t)
// where f is the firing rate function and w is the connectivity kernel
template<size_t N>
struct NeuralField {
    float tau = 10.0f;
    float threshold = 0.0f;
    float steepness = 1.0f;

    // Mexican hat connectivity: w(x) = A*exp(-x²/σ_e²) - B*exp(-x²/σ_i²)
    float A = 1.5f;      // Excitation amplitude
    float B = 0.75f;     // Inhibition amplitude
    float sigma_e = 1.0f;// Excitation width
    float sigma_i = 2.0f;// Inhibition width

    std::array<float, N> u;      // Activity field
    std::array<float, N> input;  // External input

    void initialize() {
        u.fill(0.0f);
        input.fill(0.0f);
    }

    // Firing rate function
    float f(float x) const {
        return 1.0f / (1.0f + std::exp(-steepness * (x - threshold)));
    }

    // Connectivity kernel
    float w(float dx) const {
        float dx2 = dx * dx;
        return A * std::exp(-dx2 / (2.0f * sigma_e * sigma_e)) -
               B * std::exp(-dx2 / (2.0f * sigma_i * sigma_i));
    }

    void step(float dt) {
        std::array<float, N> du;

        for (size_t i = 0; i < N; ++i) {
            // Compute integral ∫w(x-x')*f(u(x'))dx'
            float integral = 0.0f;
            for (size_t j = 0; j < N; ++j) {
                float dx = static_cast<float>(i) - static_cast<float>(j);
                integral += w(dx) * f(u[j]);
            }
            integral /= N;  // Normalize

            du[i] = (-u[i] + integral + input[i]) / tau;
        }

        for (size_t i = 0; i < N; ++i) {
            u[i] += du[i] * dt;
        }
    }
};

// Kuramoto Model (synchronization)
// dθ_i/dt = ω_i + (K/N) * Σ_j sin(θ_j - θ_i)
template<size_t N>
struct KuramotoModel {
    float K = 1.0f;  // Coupling strength

    std::array<float, N> theta;  // Phases
    std::array<float, N> omega;  // Natural frequencies

    void initialize(float omega_mean, float omega_std) {
        for (size_t i = 0; i < N; ++i) {
            // Random phase in [0, 2π]
            theta[i] = static_cast<float>(i) / N * 2.0f * M_PI;
            // Random frequency around mean
            omega[i] = omega_mean + omega_std * (static_cast<float>(i) / N - 0.5f);
        }
    }

    // Order parameter r*e^(iψ) = (1/N) * Σ_j e^(iθ_j)
    float order_parameter() const {
        float cos_sum = 0.0f, sin_sum = 0.0f;
        for (size_t i = 0; i < N; ++i) {
            cos_sum += std::cos(theta[i]);
            sin_sum += std::sin(theta[i]);
        }
        return std::sqrt(cos_sum*cos_sum + sin_sum*sin_sum) / N;
    }

    void step(float dt) {
        std::array<float, N> dtheta;

        for (size_t i = 0; i < N; ++i) {
            float coupling = 0.0f;
            for (size_t j = 0; j < N; ++j) {
                coupling += std::sin(theta[j] - theta[i]);
            }
            dtheta[i] = omega[i] + (K / N) * coupling;
        }

        for (size_t i = 0; i < N; ++i) {
            theta[i] += dtheta[i] * dt;
            // Wrap to [0, 2π]
            while (theta[i] > 2.0f * M_PI) theta[i] -= 2.0f * M_PI;
            while (theta[i] < 0.0f) theta[i] += 2.0f * M_PI;
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// LEVEL 3: BRAIN SYSTEM MATHEMATICS
// ═══════════════════════════════════════════════════════════════════════════════

// Cortical Column Model
// 6 layers: I, II/III, IV, V, VI
struct CorticalColumn {
    // Layer populations (E/I per layer)
    std::array<WilsonCowan, 6> layers;

    // Inter-layer connectivity pattern
    // Canonical microcircuit: Input → L4 → L2/3 → L5 → L6 → Thalamus
    void initialize() {
        // Set layer-specific parameters
        for (int i = 0; i < 6; ++i) {
            layers[i].E = 0.1f;
            layers[i].I = 0.1f;
        }
    }

    // Feedforward processing
    void feedforward_input(float input) {
        // Input arrives at Layer 4
        layers[3].I_E = input;
    }

    void step(float dt) {
        // L4 → L2/3
        layers[1].I_E = layers[3].E * 0.5f;

        // L2/3 → L5
        layers[4].I_E = layers[1].E * 0.6f;

        // L5 → L6
        layers[5].I_E = layers[4].E * 0.4f;

        // L6 feedback to L4
        layers[3].I_E += layers[5].E * 0.2f;

        // Update all layers
        for (auto& layer : layers) {
            layer.step(dt);
        }
    }

    // Output from L5 (to subcortical structures)
    float output() const { return layers[4].E; }
};

// Thalamo-Cortical Loop
// Thalamus acts as relay and modulator
struct ThalamoCorticalLoop {
    WilsonCowan thalamus;
    WilsonCowan cortex;
    WilsonCowan reticular;  // TRN - inhibitory

    // Connectivity
    float w_TC = 1.0f;   // Thalamus → Cortex
    float w_CT = 0.5f;   // Cortex → Thalamus
    float w_TR = 0.3f;   // Thalamus → Reticular
    float w_RT = 0.8f;   // Reticular → Thalamus (inhibitory)
    float w_CR = 0.2f;   // Cortex → Reticular

    void step(float dt, float sensory_input) {
        // Thalamus receives sensory input and cortical feedback
        thalamus.I_E = sensory_input + w_CT * cortex.E;
        thalamus.I_I = w_RT * reticular.E;  // TRN inhibition

        // Cortex receives thalamic input
        cortex.I_E = w_TC * thalamus.E;

        // Reticular receives collaterals
        reticular.I_E = w_TR * thalamus.E + w_CR * cortex.E;

        thalamus.step(dt);
        cortex.step(dt);
        reticular.step(dt);
    }
};

// Basal Ganglia Circuit
// Direct and Indirect pathways for action selection
struct BasalGanglia {
    // Nuclei
    float striatum_D1 = 0.0f;  // Direct pathway
    float striatum_D2 = 0.0f;  // Indirect pathway
    float GPe = 0.5f;          // Globus Pallidus external
    float GPi = 0.5f;          // Globus Pallidus internal
    float STN = 0.2f;          // Subthalamic Nucleus
    float SNr = 0.5f;          // Substantia Nigra reticulata

    // Dopamine modulation
    float dopamine = 1.0f;     // D1 enhance, D2 suppress

    float tau = 20.0f;

    void step(float dt, float cortical_input) {
        // Cortex → Striatum
        float d_D1 = (-striatum_D1 + cortical_input * dopamine) / tau;
        float d_D2 = (-striatum_D2 + cortical_input / dopamine) / tau;

        // Direct: Striatum D1 → GPi (inhibitory)
        // Indirect: Striatum D2 → GPe → STN → GPi
        float d_GPe = (-GPe + (1.0f - striatum_D2)) / tau;
        float d_STN = (-STN + (1.0f - GPe)) / tau;
        float d_GPi = (-GPi + STN - striatum_D1) / tau;

        striatum_D1 += d_D1 * dt;
        striatum_D2 += d_D2 * dt;
        GPe += d_GPe * dt;
        STN += d_STN * dt;
        GPi += d_GPi * dt;

        // Clamp
        striatum_D1 = std::clamp(striatum_D1, 0.0f, 1.0f);
        striatum_D2 = std::clamp(striatum_D2, 0.0f, 1.0f);
        GPe = std::clamp(GPe, 0.0f, 1.0f);
        STN = std::clamp(STN, 0.0f, 1.0f);
        GPi = std::clamp(GPi, 0.0f, 1.0f);
    }

    // Thalamic disinhibition (output)
    float output() const { return 1.0f - GPi; }
};

// Cerebellar Computation (error-driven learning)
struct Cerebellum {
    static constexpr size_t N_GRANULE = 1000;
    static constexpr size_t N_PURKINJE = 10;

    // Mossy fiber input
    std::array<float, 64> mossy_input;

    // Granule cell expansion (sparse coding)
    std::array<float, N_GRANULE> granule;

    // Purkinje cells
    std::array<float, N_PURKINJE> purkinje;

    // Parallel fiber → Purkinje weights (plastic)
    std::array<std::array<float, N_GRANULE>, N_PURKINJE> weights;

    // Climbing fiber error signal
    float climbing_fiber = 0.0f;

    float learning_rate = 0.01f;

    void initialize() {
        mossy_input.fill(0.0f);
        granule.fill(0.0f);
        purkinje.fill(0.0f);

        // Random initial weights
        for (auto& w_row : weights) {
            for (auto& w : w_row) {
                w = 0.5f;
            }
        }
    }

    void step(float dt) {
        // Granule cell activation (sparse)
        for (size_t i = 0; i < N_GRANULE; ++i) {
            float input_sum = 0.0f;
            for (size_t j = 0; j < mossy_input.size(); ++j) {
                input_sum += mossy_input[j] * ((i * j) % 2 == 0 ? 1.0f : 0.0f);
            }
            granule[i] = (input_sum > 0.5f) ? 1.0f : 0.0f;
        }

        // Purkinje cell computation
        for (size_t p = 0; p < N_PURKINJE; ++p) {
            float pf_input = 0.0f;
            for (size_t g = 0; g < N_GRANULE; ++g) {
                pf_input += weights[p][g] * granule[g];
            }
            purkinje[p] = 1.0f / (1.0f + std::exp(-pf_input + 5.0f));
        }

        // LTD at parallel fiber synapses (climbing fiber signal)
        if (climbing_fiber > 0.5f) {
            for (size_t p = 0; p < N_PURKINJE; ++p) {
                for (size_t g = 0; g < N_GRANULE; ++g) {
                    if (granule[g] > 0.5f) {
                        weights[p][g] -= learning_rate * climbing_fiber;
                        weights[p][g] = std::max(0.0f, weights[p][g]);
                    }
                }
            }
        }
    }

    // Deep cerebellar nuclei output (disinhibition)
    float output() const {
        float sum = 0.0f;
        for (float p : purkinje) sum += p;
        return 1.0f - sum / N_PURKINJE;
    }
};

// Spinal Cord Reflex Arc
// Simple 3-neuron circuit: Sensory → Interneuron → Motor
struct SpinalReflex {
    HodgkinHuxley sensory;
    HodgkinHuxley interneuron;
    HodgkinHuxley motor;

    AMPASynapse s_to_i;  // Sensory → Interneuron
    AMPASynapse i_to_m;  // Interneuron → Motor

    float sensory_threshold = 5.0f;

    void step(float dt, float stimulus) {
        // Sensory neuron responds to stimulus
        sensory.step(dt, stimulus);

        // If sensory spikes, activate synapse
        if (sensory.is_spiking()) {
            s_to_i.receive_spike();
        }
        s_to_i.step(dt);

        // Interneuron receives synaptic current
        interneuron.step(dt, s_to_i.current(interneuron.V));

        // Interneuron drives motor neuron
        if (interneuron.is_spiking()) {
            i_to_m.receive_spike();
        }
        i_to_m.step(dt);

        motor.step(dt, i_to_m.current(motor.V));
    }

    // Motor output (muscle activation)
    float motor_output() const {
        return motor.is_spiking() ? 1.0f : 0.0f;
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// INTEGRATED NERVOUS SYSTEM
// Connects all components
// ═══════════════════════════════════════════════════════════════════════════════

struct NervousSystem {
    // Cortical areas
    std::array<CorticalColumn, 4> cortical_columns;  // Simplified cortex

    // Subcortical structures
    ThalamoCorticalLoop thalamo_cortical;
    BasalGanglia basal_ganglia;
    Cerebellum cerebellum;

    // Spinal cord
    std::array<SpinalReflex, 4> spinal_segments;

    // Global state
    float arousal = 0.5f;      // Reticular activating system
    float attention = 0.5f;    // Selective attention modulation

    void initialize() {
        for (auto& col : cortical_columns) col.initialize();
        cerebellum.initialize();
    }

    void step(float dt, const std::array<float, 4>& sensory_inputs) {
        // 1. Sensory input → Thalamus → Cortex
        float avg_input = 0.0f;
        for (float s : sensory_inputs) avg_input += s;
        avg_input /= 4.0f;

        thalamo_cortical.step(dt, avg_input * arousal);

        // 2. Cortex → Basal Ganglia (action selection)
        basal_ganglia.step(dt, thalamo_cortical.cortex.E);

        // 3. Cortex + Basal Ganglia → Cerebellum (motor learning)
        for (size_t i = 0; i < cerebellum.mossy_input.size(); ++i) {
            cerebellum.mossy_input[i] = (i < 4) ? sensory_inputs[i] : 0.0f;
        }
        cerebellum.step(dt);

        // 4. Motor output → Spinal cord
        float motor_command = basal_ganglia.output() * cerebellum.output();
        for (size_t i = 0; i < 4; ++i) {
            // Descending motor command + local sensory
            float stimulus = motor_command * 10.0f + sensory_inputs[i] * 5.0f;
            spinal_segments[i].step(dt, stimulus);
        }

        // 5. Update arousal (homeostatic)
        arousal = 0.99f * arousal + 0.01f * avg_input;
    }

    // Get motor outputs
    std::array<float, 4> motor_outputs() const {
        std::array<float, 4> out;
        for (size_t i = 0; i < 4; ++i) {
            out[i] = spinal_segments[i].motor_output();
        }
        return out;
    }

    // Consciousness metric (integrated information)
    float consciousness() const {
        float phi = 0.0f;
        phi += thalamo_cortical.cortex.E;
        phi += thalamo_cortical.thalamus.E;
        phi += (1.0f - basal_ganglia.GPi);
        phi += arousal;
        return phi / 4.0f;
    }
};

} // namespace neuro
} // namespace mind_farm