% =========================================================================
% Optimising propofol induction in elderly patients
% by constraining peak plasma concentration or infusion rate during 
% Eleveld target-controlled infusions: an in silico simulation study
%
% George Zhong @ Propofol Dreams
%
% Compares two target-controlled infusion (TCI) induction strategies to standard
% effect-site targeting:
%   - Set A (Cp-Limited): Infuses at max pump rate until target peak plasma
%     concentration (Cp_lim) is reached, maintain Cp for time t_clamp, then coasts.
%   - Set B (Rate-Limited): Constant infusion at rate 'p' until pause,
%     then coasts until reaching target Ce.
%
% -------------------------------------------------------------------------
% ALGORITHM OVERVIEW
% -------------------------------------------------------------------------
% PK model (Eleveld 3-compartment, amounts + effect site):
%   State x = [A1; A2; A3; Ce]  (A1..A3 in mg, Ce in mcg/mL)
%   A_amt is a 4x4 linear time-invariant system matrix; B_amt selects the
%   infusion input into the central compartment.
%
% Analytical solvers (see functions at bottom of file):
%   compute_state_free(t,x0)   -> expm(A*t) * x0          (coasting)
%   compute_state_forced(t,u)  -> A\(expm(A*t)-I)*B*u     (constant infusion)
%   find_free_peak_Ce(x0)      -> solves Cp(t) = Ce(t)    (Ce peak on coast)
%   get_state_after_clamp_fast -> closed-form peripheral/effect-site
%                                 redistribution while Cp is clamped.
%
% Set A (Cp-Limited) phases:
%   Phase 1 (ramp):   u = max_rate until Cp reaches Cp_lim.
%   Phase 2 (clamp):  u chosen algebraically so dA1/dt = 0, i.e.
%                     u = k_elim*Cp_lim*V1 - k21*A2 - k31*A3,
%                     keeping A1 = Cp_lim*V1. Peripheral/effect-site
%                     states relax exponentially toward their equilibria.
%   Phase 3 (coast):  u = 0 until Ce peaks at target_Ce.
%   Phase 4 (TCI):    u >= 0 chosen to maintain Ce at target_Ce thereafter.
%
% Set B (Rate-Limited) phases:
%   Phase 1 (infusion): u = p (constant) until Ce peaks at target_Ce.
%   Phase 2 (TCI):      u >= 0 chosen to maintain Ce = target_Ce.
%
% Feasibility sweeps map the parameter 'p' (Cp_lim for A, infusion rate for
% B) to the resulting time-to-meet via a monotonic log-spaced grid. The
% feasible window [t_min, t_max] is the set of achievable target meet times.
% =========================================================================

clear; clc; close all;

%% 1. PATIENT COVARIATES & ELEVELD PK-PD PARAMETERIZATION
% Patient Profile: 80yo Male, 70kg, 170cm
age    = 80;    % years
weight = 70;    % kg
height = 170;   % cm
sex    = 1;     % 1 = Male, 2 = Female

% Central and shallow peripheral distribution volumes (L)
V1 = (6.28 * (weight / (weight + 33.6))) / 0.675675675676;
V2 = 25.5 * (weight / 70) * exp(-0.0156 * (age - 35));

hMeters  = height / 100;
bmi_term = weight / (hMeters^2);

% Lean Body Mass (LBM) using Janmahasatian formulations in Eleveld model
if sex == 1
    patientLBM = (0.88 + 0.12 / (1 + (age / 13.4)^-12.7)) * ((9270 * weight) / (6680 + 216 * bmi_term));
    cl1_ref    = 1.79;
else
    patientLBM = (1.11 - 0.11 / (1 + (age / 7.1)^-1.1)) * ((9270 * weight) / (8780 + 244 * bmi_term));
    cl1_ref    = 2.10;
end

% Deep peripheral volume (L) and Post-Menstrual Age (PMA) maturation
V3      = (273 * exp(-0.0138 * age) * patientLBM) / 54.4752059601377;
pma     = age * 52.143 + 40;
pma_mat = (pma^9.06) / (pma^9.06 + 42.3^9.06);

% Clearances (L/min)
Cl1     = cl1_ref * (weight / 70)^0.75 * pma_mat * exp(-0.00286 * age);
v2_term = (25.5 * (weight / 70) * exp(-0.0156 * (age - 35))) / 25.5;
Cl2     = 1.75 * v2_term^0.75 * (1 + 1.3 * (1 - pma / (pma + 68.3)));
v3_norm = (patientLBM * exp(-0.0138 * age)) / 54.4752059601377;
Cl3     = 1.11 * v3_norm^0.75 * ((pma / (pma + 68.3)) / 0.964695544);

% Micro-rate constants (converted from 1/min to 1/s)
k10 = (Cl1 / V1) / 60;
k12 = (Cl2 / V1) / 60;
k21 = (Cl2 / V2) / 60;
k13 = (Cl3 / V1) / 60;
k31 = (Cl3 / V3) / 60;
ke0 = (0.146 * (weight / 70)^-0.25) / 60; % Effect-site transfer rate constant

% Infusion Limits & Targets
max_rate  = (1200 * 10) / 3600; % Max pump speed: 1200 mL/h @ 10 mg/mL = 3.33 mg/s
target_Ce = 3.0;                % Target effect-site concentration (mcg/mL)
dt        = 0.1;                % Simulation time step (s)

%% 2. STATE-SPACE MATRICES & SPECTRAL ACCELERATION
% State vector: x = [A1 (mg); A2 (mg); A3 (mg); Ce (mcg/mL)]
A_amt = [ -(k10+k12+k13),  k21,  k31,  0;
           k12,           -k21,  0,    0;
           k13,            0,   -k31,  0;
           ke0/V1,         0,    0,   -ke0 ];
B_amt = [1; 0; 0; 0];

% Precompute eigendecomposition: expm(A*t) = V_eig * diag(exp(lambda*t)) * invV
[V_eig, D_eig] = eig(A_amt);
lambda         = diag(D_eig);
invV           = inv(V_eig);

% Package model system cache for fast evaluations
sys = struct('A', A_amt, 'B', B_amt, 'V', V_eig, 'lambda', lambda, ...
             'invV', invV, 'V1', V1, 'ke0', ke0, 'k10', k10, ...
             'k12', k12, 'k13', k13, 'k21', k21, 'k31', k31, ...
             'target_Ce', target_Ce, 'max_rate', max_rate);

%% 3. FEASIBILITY SWEEP: BOUNDARY IDENTIFICATION
disp('Calibrating physical boundaries (dt = 0.1s)...');

% --- Set B: Rate-limited boundary sweep ---
steady_state_elim = k10 * target_Ce * V1;
p_B_lo   = steady_state_elim * 1.0001;
p_B_hi   = max_rate;
n_sweep  = 200;

p_B_sweep = logspace(log10(p_B_lo), log10(p_B_hi), n_sweep);
t_B_sweep = nan(size(p_B_sweep));
for i = 1:n_sweep
    t_B_sweep(i) = get_t_meet_fast(2, p_B_sweep(i), sys);
end

feas_B   = isfinite(t_B_sweep) & t_B_sweep > 0 & t_B_sweep < 3600;
p_B_feas = p_B_sweep(feas_B);
t_B_feas = t_B_sweep(feas_B);
if isempty(t_B_feas)
    error('Set B (rate-limited) cannot reach target_Ce for this patient.');
end
t_B_min = min(t_B_feas);
t_B_max = max(t_B_feas);

% --- Set A: Cp-limited boundary sweep ---
obj_inf_max = @(ti) eval_peak_Ce_inf_fast(ti, max_rate, sys) - target_Ce;
if obj_inf_max(0) >= 0
    t_inf_max = 0;
elseif obj_inf_max(3600) < 0
    error('max_rate cannot reach target_Ce within 3600s.');
else
    t_inf_max = fzero(obj_inf_max, [0, 3600]);
end

x_max_inf = compute_state_forced(t_inf_max, max_rate, sys);
Max_CpT   = x_max_inf(1) / V1;
max_y_cp  = ceil(Max_CpT * 1.1);

p_A_lo    = target_Ce * 1.001;
p_A_hi    = Max_CpT;
p_A_sweep = logspace(log10(p_A_lo), log10(p_A_hi), n_sweep);
t_A_sweep = nan(size(p_A_sweep));
for i = 1:n_sweep
    t_A_sweep(i) = get_t_meet_fast(1, p_A_sweep(i), sys);
end

feas_A   = isfinite(t_A_sweep) & t_A_sweep > 0 & t_A_sweep < 3600;
p_A_feas = p_A_sweep(feas_A);
t_A_feas = t_A_sweep(feas_A);
if isempty(t_A_feas)
    error('Set A (Cp-limited) cannot reach target_Ce for this patient.');
end
t_A_min = min(t_A_feas);
t_A_max = max(t_A_feas);

% Overlapping feasible induction window
t_common_min    = max(t_A_min, t_B_min);
target_max_time = 900;
t_common_max    = min(min(t_A_max, t_B_max), target_max_time);

if t_common_min >= t_common_max
    error('No overlapping feasible target-time window between strategies within 0-900s.');
end

fprintf('Set A feasible t_meet range: [%.1f, %.1f] s\n', t_A_min, t_A_max);
fprintf('Set B feasible t_meet range: [%.1f, %.1f] s\n', t_B_min, t_B_max);
fprintf('Common feasible range (capped at 900s): [%.1f, %.1f] s\n', t_common_min, t_common_max);

%% 4. BASELINE REFERENCE & RUN CONFIGURATION
[baseline_t_meet, base_Cp, ~, base_Rate] = simulate_regime_fast(2, max_rate, sys, dt);
baseline_peak_cp = max(base_Cp);
end_idx_base     = min(round(900 / dt), length(base_Rate));
baseline_total_dose = sum(base_Rate(1:end_idx_base)) * dt;

% Discretized target times: 20 points from fastest meet time to 900s
n_targets       = 20;
planned_targets = linspace(t_common_min, target_max_time, n_targets);
time_vector     = dt:dt:1800;
n_sim_pts       = length(time_vector);

% Preallocate result structures
res = struct();
for s = 1:2
    res(s).t_meet        = zeros(1, n_targets);
    res(s).peak_cp       = zeros(1, n_targets);
    res(s).dose_15       = zeros(1, n_targets);
    res(s).Cp_profiles   = zeros(n_targets, n_sim_pts);
    res(s).Ce_profiles   = zeros(n_targets, n_sim_pts);
    res(s).Rate_profiles = zeros(n_targets, n_sim_pts);
end

%% 5. EXECUTION & SIMULATION SWEEP
disp('------------------------------------------------------');
disp(' Target Meet | Set A Actual (Diff) | Set B Actual (Diff)');
disp('------------------------------------------------------');

for s = 1:2
    if s == 1
        p_lo = p_A_lo; p_hi = p_A_hi;
    else
        p_lo = p_B_lo; p_hi = p_B_hi;
    end

    for idx = 1:n_targets
        tt = planned_targets(idx);
        if idx == 1
            best_mid = p_hi;
        else
            best_mid = solve_bracketed_zero(@(p) get_t_meet_fast(s, p, sys) - tt, p_lo, p_hi);
        end
        [final_t, Cp_arr, Ce_arr, Rate_arr] = simulate_regime_fast(s, best_mid, sys, dt);

        res(s).t_meet(idx)          = final_t;
        res(s).peak_cp(idx)         = max(Cp_arr);
        end_idx                     = min(round(900 / dt), n_sim_pts);
        res(s).dose_15(idx)         = sum(Rate_arr(1:end_idx)) * dt;
        res(s).Cp_profiles(idx,:)   = Cp_arr;
        res(s).Ce_profiles(idx,:)   = Ce_arr;
        res(s).Rate_profiles(idx,:) = Rate_arr;
    end
end

% Display tabular differences
for i = 1:n_targets
    diff_A = abs(planned_targets(i) - res(1).t_meet(i));
    diff_B = abs(planned_targets(i) - res(2).t_meet(i));
    fprintf(' %7.1fs   |     %7.1fs (%4.1fs)  |     %7.1fs (%4.1fs)\n', ...
        planned_targets(i), res(1).t_meet(i), diff_A, res(2).t_meet(i), diff_B);
end
disp('------------------------------------------------------');

%% 6. VISUALIZATION & EXPORT
% Set explicit physical dimensions in inches so -r600 produces exact pixel dimensions
% and embeds strict 600 DPI metadata in the PNG/TIFF header.
fig1 = figure('Color', 'w', 'Name', 'Concentration & Infusion Profiles', ...
              'Units', 'inches', 'Position', [0.5, 0.5, 12, 10], ...
              'PaperUnits', 'inches', 'PaperPosition', [0, 0, 12, 10], ...
              'PaperSize', [12, 10]);

fig2 = figure('Color', 'w', 'Name', 'Induction Analysis', ...
              'Units', 'inches', 'Position', [1.0, 1.0, 12, 5], ...
              'PaperUnits', 'inches', 'PaperPosition', [0, 0, 12, 5], ...
              'PaperSize', [12, 5]);

titles       = {'Peak Cp limited', 'Infusion rate limited'};
group_colors = [0.5 0.5 0.5; 0 0 0];

% Figure 1: Concentration & Rate Over Time (Panels A, B, C, D)
panel_labels_fig1 = {'A', 'B', 'C', 'D'};
figure(fig1);
for s = 1:2
    % Concentration subplots (Cp and Ce)
    subplot(2, 2, s); hold on;
    plot(time_vector, res(s).Cp_profiles, 'Color', [0.8 0.8 0.8], 'LineWidth', 1.0);
    plot(time_vector, res(s).Ce_profiles, 'Color', [0.8 0.8 0.8], 'LineWidth', 1.5, 'LineStyle', '--');
    plot(time_vector, res(s).Cp_profiles(3,:), 'Color', 'k', 'LineWidth', 1.5);
    plot(time_vector, res(s).Ce_profiles(3,:), 'Color', 'k', 'LineWidth', 2.0, 'LineStyle', '--');
    xlabel('Time (s)'); ylabel('Concentration (\mug\cdotmL^{-1})');
    xlim([0 900]); xticks(0:100:900);
    apply_axis_style(gca);
    text(-0.12, 1.06, panel_labels_fig1{s}, 'Units', 'normalized', ...
         'FontSize', 14, 'FontWeight', 'bold', 'Color', 'k');

    % Infusion Rate subplots
    subplot(2, 2, s+2); hold on;
    plot(time_vector, res(s).Rate_profiles .* 360, 'Color', [0.8 0.8 0.8], 'LineWidth', 1.5);
    plot(time_vector, res(s).Rate_profiles(3,:) .* 360, 'Color', 'k', 'LineWidth', 2.0);
    xlabel('Time (s)'); ylabel('Pump Rate (mL\cdoth^{-1})');
    xlim([0 900]); xticks(0:100:900);
    apply_axis_style(gca);
    text(-0.12, 1.06, panel_labels_fig1{s+2}, 'Units', 'normalized', ...
         'FontSize', 14, 'FontWeight', 'bold', 'Color', 'k');
end

% Figure 2: Speed vs. Haemodynamic Surrogates (Panels A and B)
panel_labels_fig2 = {'A', 'B'};
figure(fig2);
for s = 1:2
    subplot(1, 2, 1); hold on;
    y1_plot = (res(s).peak_cp ./ baseline_peak_cp) .* 100;
    scatter(res(s).t_meet, y1_plot, 80, group_colors(s,:), 'filled', ...
            'MarkerEdgeColor', 'k', 'DisplayName', titles{s});

    subplot(1, 2, 2); hold on;
    y2_plot = (res(s).dose_15 ./ baseline_total_dose) .* 100;
    scatter(res(s).t_meet, y2_plot, 80, group_colors(s,:), 'filled', ...
            'MarkerEdgeColor', 'k', 'DisplayName', titles{s});
end

subplot(1, 2, 1);
xlabel('Time to Ce target (s)'); ylabel('Peak Cp (% of baseline)');
title('Induction Speed vs. Peak Cp');
lgd1 = legend(titles, 'Location', 'best');
set(lgd1, 'Color', 'w', 'EdgeColor', 'k', 'TextColor', 'k', 'FontSize', 10);
apply_axis_style(gca);
text(-0.12, 1.06, panel_labels_fig2{1}, 'Units', 'normalized', ...
     'FontSize', 14, 'FontWeight', 'bold', 'Color', 'k');

subplot(1, 2, 2);
xlabel('Time to Ce target (s)'); ylabel('15 Min Dose (% of baseline)');
title('Induction Speed vs. 15-Min Propofol Dose');
lgd2 = legend(titles, 'Location', 'best');
set(lgd2, 'Color', 'w', 'EdgeColor', 'k', 'TextColor', 'k', 'FontSize', 10);
apply_axis_style(gca);
text(-0.12, 1.06, panel_labels_fig2{2}, 'Units', 'normalized', ...
     'FontSize', 14, 'FontWeight', 'bold', 'Color', 'k');

% Export Publication Graphics at strict 600 DPI
drawnow;
print(fig1, 'Figure1.png', '-dpng', '-r600');
print(fig2, 'Figure2.png', '-dpng', '-r600');

% Optional: 600 DPI TIFF export (uncomment if journal requires TIFF)
% print(fig1, 'Figure1.tif', '-dtiff', '-r600');
% print(fig2, 'Figure2.tif', '-dtiff', '-r600');

%% ========================================================================
%  ANALYTICAL PHARMACOKINETIC SOLVER FUNCTIONS
%  ========================================================================

function [t_meet, Cp_arr, Ce_arr, Rate_arr] = simulate_regime_fast(set_id, p, sys, dt)
    t_meet = get_t_meet_fast(set_id, p, sys);

    if set_id == 1
        Cp_lim = p;
        obj_up = @(t) ([1/sys.V1, 0, 0, 0] * compute_state_forced(t, sys.max_rate, sys)) - Cp_lim;
        t_up   = solve_bracketed_zero(obj_up, 0, 3600);
        x_clamp_start = compute_state_forced(t_up, sys.max_rate, sys);

        obj_cl = @(tc) eval_peak_Ce_clamp_fast(tc, Cp_lim, x_clamp_start, sys) - sys.target_Ce;
        t_clamp = solve_bracketed_zero(obj_cl, 0, 3600);
        t_inf   = 0;
    else
        Cp_lim  = 0;
        t_up    = 0;
        t_clamp = 0;
        obj_inf = @(ti) eval_peak_Ce_inf_fast(ti, p, sys) - sys.target_Ce;
        t_inf   = solve_bracketed_zero(obj_inf, 0, 3600);
    end

    steps    = round(1800 / dt);
    Cp_arr   = zeros(1, steps);
    Ce_arr   = zeros(1, steps);
    Rate_arr = zeros(1, steps);

    Phi   = sys.V * diag(exp(sys.lambda * dt)) * sys.invV;
    Gamma = (sys.A \ (Phi - eye(4))) * sys.B;
    x     = zeros(4, 1);

    exp_k21 = exp(-sys.k21 * dt);
    exp_k31 = exp(-sys.k31 * dt);
    exp_ke0 = exp(-sys.ke0 * dt);

    k_elim = sys.k10 + sys.k12 + sys.k13;

    Cp_lim_V1 = Cp_lim * sys.V1;
    X2_eq_cl  = (sys.k12 / sys.k21) * Cp_lim_V1;
    X3_eq_cl  = (sys.k13 / sys.k31) * Cp_lim_V1;
    Tc_V1     = sys.target_Ce * sys.V1;
    X2_eq_tg  = (sys.k12 / sys.k21) * Tc_V1;
    X3_eq_tg  = (sys.k13 / sys.k31) * Tc_V1;

    clamp_end = t_up + t_clamp;

    for step = 1:steps
        t = step * dt;

        if set_id == 1
            if t <= t_up
                u = sys.max_rate;
            elseif t <= clamp_end
                u = max(0, k_elim * Cp_lim_V1 - sys.k21 * x(2) - sys.k31 * x(3));
            elseif t <= t_meet
                u = 0;
            else
                u = max(0, k_elim * Tc_V1 - sys.k21 * x(2) - sys.k31 * x(3));
            end
        else
            if t <= t_inf
                u = p;
            elseif t <= t_meet
                u = 0;
            else
                u = max(0, k_elim * Tc_V1 - sys.k21 * x(2) - sys.k31 * x(3));
            end
        end
        Rate_arr(step) = u;

        if set_id == 1 && t > t_up && t <= clamp_end
            x(1) = Cp_lim_V1;
            x(2) = X2_eq_cl + (x(2) - X2_eq_cl) * exp_k21;
            x(3) = X3_eq_cl + (x(3) - X3_eq_cl) * exp_k31;
            x(4) = Cp_lim    + (x(4) - Cp_lim)    * exp_ke0;
        elseif t > t_meet
            x(1) = Tc_V1;
            x(2) = X2_eq_tg + (x(2) - X2_eq_tg) * exp_k21;
            x(3) = X3_eq_tg + (x(3) - X3_eq_tg) * exp_k31;
            x(4) = sys.target_Ce + (x(4) - sys.target_Ce) * exp_ke0;
        else
            x = Phi * x + Gamma * u;
        end

        Cp_arr(step) = x(1) / sys.V1;
        Ce_arr(step) = x(4);
    end
end

function t_meet = get_t_meet_fast(set_id, p, sys)
    if set_id == 1
        Cp_lim = p;
        obj_up = @(t) ([1/sys.V1, 0, 0, 0] * compute_state_forced(t, sys.max_rate, sys)) - Cp_lim;
        t_up   = solve_bracketed_zero(obj_up, 0, 3600);
        if isnan(t_up), t_meet = NaN; return; end

        x_clamp_start = compute_state_forced(t_up, sys.max_rate, sys);

        obj_cl = @(tc) eval_peak_Ce_clamp_fast(tc, Cp_lim, x_clamp_start, sys) - sys.target_Ce;
        t_clamp = solve_bracketed_zero(obj_cl, 0, 3600);
        if isnan(t_clamp), t_meet = NaN; return; end

        x_coast_start = get_state_after_clamp_fast(t_clamp, Cp_lim, x_clamp_start, sys);
        obj_cst = @(tcst) ([1/sys.V1, 0, 0, 0] * compute_state_free(tcst, x_coast_start, sys)) - sys.target_Ce;
        t_coast = solve_bracketed_zero(obj_cst, 0, 3600);
        if isnan(t_coast), t_meet = NaN; return; end

        t_meet = t_up + t_clamp + t_coast;
    else
        obj_inf = @(ti) eval_peak_Ce_inf_fast(ti, p, sys) - sys.target_Ce;
        t_inf   = solve_bracketed_zero(obj_inf, 0, 3600);
        if isnan(t_inf), t_meet = NaN; return; end

        x_coast_start = compute_state_forced(t_inf, p, sys);
        obj_cst = @(tcst) ([1/sys.V1, 0, 0, 0] * compute_state_free(tcst, x_coast_start, sys)) - sys.target_Ce;
        t_coast = solve_bracketed_zero(obj_cst, 0, 3600);
        if isnan(t_coast), t_meet = NaN; return; end

        t_meet = t_inf + t_coast;
    end
end

%% --- EXACT PK ANALYTICAL EVALUATORS ---

function x_t = compute_state_free(t, x0, sys)
    x_t = sys.V * (exp(sys.lambda * t) .* (sys.invV * x0));
end

function x_t = compute_state_forced(t, rate, sys)
    exp_term = exp(sys.lambda * t) - 1;
    invV_B   = sys.invV * (sys.B * rate);
    int_exp  = sys.V * (exp_term .* (invV_B ./ sys.lambda));
    x_t      = int_exp;
end

function peak = eval_peak_Ce_clamp_fast(tc, Cp_lim, x_start, sys)
    x_coast = get_state_after_clamp_fast(tc, Cp_lim, x_start, sys);
    peak    = find_free_peak_Ce(x_coast, sys);
end

function peak = eval_peak_Ce_inf_fast(t_inf, rate, sys)
    x_coast = compute_state_forced(t_inf, rate, sys);
    peak    = find_free_peak_Ce(x_coast, sys);
end

function peak = find_free_peak_Ce(x0, sys)
    dCe_dt = @(t) ([1/sys.V1, 0, 0, -1] * compute_state_free(t, x0, sys));

    if dCe_dt(0) <= 0
        t_peak = 0;
    else
        t_peak = solve_bracketed_zero(dCe_dt, 0, 1800);
        if isnan(t_peak), t_peak = 0; end
    end
    x_peak = compute_state_free(t_peak, x0, sys);
    peak   = x_peak(4);
end

function x = get_state_after_clamp_fast(t_clamp, Cp_lim, x_start, sys)
    X1    = Cp_lim * sys.V1;
    X2_eq = (sys.k12 / sys.k21) * X1;
    X3_eq = (sys.k13 / sys.k31) * X1;

    x    = zeros(4, 1);
    x(1) = X1;
    x(2) = X2_eq + (x_start(2) - X2_eq) * exp(-sys.k21 * t_clamp);
    x(3) = X3_eq + (x_start(3) - X3_eq) * exp(-sys.k31 * t_clamp);
    x(4) = Cp_lim + (x_start(4) - Cp_lim) * exp(-sys.ke0 * t_clamp);
end

function root = solve_bracketed_zero(fun, lower_b, upper_b)
    f_low = fun(lower_b);
    if f_low >= 0 && fun(upper_b) >= 0
        root = lower_b; return;
    elseif f_low <= 0 && fun(upper_b) <= 0
        root = upper_b; return;
    end
    try
        root = fzero(fun, [lower_b, upper_b]);
    catch
        root = NaN;
    end
end

function apply_axis_style(ax)
    set(ax, 'Box', 'off', 'TickDir', 'out', 'LineWidth', 1.5, ...
            'FontSize', 12, 'Color', 'w', 'XColor', 'k', 'YColor', 'k', ...
            'FontWeight', 'bold');

    if ~isempty(ax.Title) && isprop(ax.Title, 'Color')
        ax.Title.Color = 'k';
        ax.Title.FontWeight = 'bold';
    end

    if ~isempty(ax.XLabel), ax.XLabel.Color = 'k'; end
    if ~isempty(ax.YLabel), ax.YLabel.Color = 'k'; end
end