DRY LAB : CODING

import numpy as np
from scipy import stats

def hill_equation(x, K, n):
    \"\"\"Hill响应函数\"\"\"
    return (x**n) / (K**n + x**n)

def and_gate_score(x1, x2, x3):
    \"\"\"三钥匙AND门得分\"\"\"
    f1 = hill_equation(x1, K1=30, n1=1.5)  # 胆汁酸
    f2 = hill_equation(x2, K2=2, n2=2.0)   # FFA
    f3 = hill_equation(x3, K3=50, n3=1.8)  # TNF-α
    # 加权求和
    S = 0.2*f1 + 0.5*f2 + 1.0*f3
    return S

# 健康人分布参数
healthy_params = {
    'bile_acid': {'mean_log': 2.5, 'sigma_log': 0.3},
    'ffa': {'mean_log': -0.5, 'sigma_log': 0.4},
    'tnf': {'mean_log': 2.8, 'sigma_log': 0.5}
}

# NASH患者分布参数
nash_params = {
    'bile_acid': {'mean_log': 3.8, 'sigma_log': 0.4},
    'ffa': {'mean_log': 1.4, 'sigma_log': 0.5},
    'tnf': {'mean_log': 4.5, 'sigma_log': 0.6}
}

def generate_samples(params, n=100000):
    \"\"\"生成对数正态分布样本\"\"\"
    bile = np.random.lognormal(params['bile_acid']['mean_log'],
                                params['bile_acid']['sigma_log'], n)
    ffa = np.random.lognormal(params['ffa']['mean_log'],
                               params['ffa']['sigma_log'], n)
    tnf = np.random.lognormal(params['tnf']['mean_log'],
                               params['tnf']['sigma_log'], n)
    return bile, ffa, tnf

# 生成样本
np.random.seed(42)
h_bile, h_ffa, h_tnf = generate_samples(healthy_params)
n_bile, n_ffa, n_tnf = generate_samples(nash_params)

# 计算AND门得分
h_scores = and_gate_score(h_bile, h_ffa, h_tnf)
n_scores = and_gate_score(n_bile, n_ffa, n_tnf)

# 激活判断（阈值=1.0）
threshold = 1.0
h_activate = np.mean(h_scores > threshold)
n_activate = np.mean(n_scores > threshold)

print(f"健康人误激活率: {h_activate:.4f} ({h_activate*100:.2f}%)")
print(f"NASH患者激活率: {n_activate:.4f} ({n_activate*100:.2f}%)")

# ROC曲线计算
from sklearn.metrics import roc_curve, auc

y_true = np.concatenate([np.zeros(len(h_scores)), np.ones(len(n_scores))])
y_scores = np.concatenate([h_scores, n_scores])
fpr, tpr, thresholds = roc_curve(y_true, y_scores)
roc_auc = auc(fpr, tpr)

print(f"\nAUC = {roc_auc:.4f}")

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def dual_drug_model(y, t, params):
    \"\"\"
    双药动力学模型
    y = [N, A, R, T, F]
    \"\"\"
    N, A, R, T, F = y
    mu, Nmax, alpha_A, beta_A, gamma_A, n, KA, alpha_R, beta_R = params[:9]
    k_Tsyn, k_Tdeg, k_Fsyn, k_Fdeg, w_T, w_F = params[9:]

    # 菌生长（Logistic）
    dN = mu * N * (1 - N/Nmax) - 0.01 * N  # 自然死亡

    # AHL合成与降解
    LuxI = 1.0  # 假设LuxI表达恒定
    dA = alpha_A * (N/1e8) * LuxI - beta_A * A - gamma_A * A

    # LuxR-AHL复合物
    dR = alpha_R * (A**n) / (KA**n + A**n) - beta_R * R

    # 资源分配权重
    total_weight = w_T + w_F
    frac_T = w_T / total_weight
    frac_F = w_F / total_weight

    # T2合成（受R激活）
    dT = k_Tsyn * R * frac_T - k_Tdeg * T
    # FGF21合成
    dF = k_Fsyn * R * frac_F - k_Fdeg * F

    return [dN, dA, dR, dT, dF]

# 参数
params = [
    0.35,        # mu
    1e9,         # Nmax
    5.0,         # alpha_A
    0.1,         # beta_A
    0.05,        # gamma_A
    2.0,         # n (Hill)
    50.0,        # KA
    10.0,        # alpha_R
    0.5,         # beta_R
    0.8,         # k_Tsyn
    0.08,        # k_Tdeg
    0.4,         # k_Fsyn
    0.08,        # k_Fdeg
    0.8,         # w_T
    0.4          # w_F
]

# 初始条件
y0 = [1e6, 0, 0, 0, 0]  # 初始菌密度1e6，其他为0
t = np.linspace(0, 72, 1000)  # 72小时

# 求解
sol = odeint(dual_drug_model, y0, t, args=(params,))
N, A, R, T, F = sol.T

# 绘图
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 菌密度
ax1 = axes[0, 0]
ax1.semilogy(t, N, 'b-', linewidth=2)
ax1.axhline(y=1e8, color='r', linestyle='--', label='Steady state (10^8)')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel('Cell density (CFU/mL)')
ax1.set_title('Cell Growth Dynamics')
ax1.legend()

# AHL和LuxR
ax2 = axes[0, 1]
ax2.plot(t, A, 'g-', label='AHL', linewidth=2)
ax2.plot(t, R*50, 'm--', label='LuxR (×50)', linewidth=2)
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Concentration (nM)')
ax2.set_title('Quorum Sensing Activation')
ax2.legend()

# 双药产量
ax3 = axes[1, 0]
ax3.plot(t, T, 'b-', linewidth=2, label=f'T2 (target={T[-1]:.1f} mg/L)')
ax3.plot(t, F, 'r-', linewidth=2, label=f'FGF21 (target={F[-1]:.1f} mg/L)')
ax3.axhline(y=10, color='b', linestyle='--', alpha=0.5)
ax3.axhline(y=5, color='r', linestyle='--', alpha=0.5)
ax3.set_xlabel('Time (h)')
ax3.set_ylabel('Concentration (mg/L)')
ax3.set_title('Dual Drug Production')
ax3.legend()

# 比例验证
ax4 = axes[1, 1]
ratio = T / (F + 1e-10)  # 避免除零
ax4.plot(t, ratio, 'k-', linewidth=2)
ax4.axhline(y=2, color='r', linestyle='--', label='Target ratio 2:1')
ax4.set_xlabel('Time (h)')
ax4.set_ylabel('T2/FGF21 ratio')
ax4.set_title('Production Ratio Stability')
ax4.set_ylim(0, 5)
ax4.legend()

plt.tight_layout()
plt.savefig('dual_drug_kinetics_math.png', dpi=300)

print(f"稳态T2: {T[-1]:.2f} mg/L")
print(f"稳态FGF21: {F[-1]:.2f} mg/L")
print(f"实际比例: {T[-1]/F[-1]:.2f}:1")
print(f"达到90%稳态时间: {t[np.where(T > 0.9*T[-1])[0][0]]:.1f} h")

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 单锁失效概率（保守估计）
p_locks = np.array([1e-3, 1e-4, 1e-5, 1e-6])
lock_names = ['ΔdapA\n(Auxotroph)', 'TA System\n(Generation Limit)', 
              'pH Lysis\n(Environmental)', 'Horizontal Blocking\n(Gene Transfer)']

# 理论联合失效概率（独立假设）
p_theory = np.prod(p_locks)
print(f"理论联合失效概率: {p_theory:.2e}")

# ==================== Monte Carlo模拟 ====================

def monte_carlo_safety(p_locks, n_simulations, n_repeats=10):
    """
    Monte Carlo估计联合失效概率
    """
    estimates = []
    
    for _ in range(n_repeats):
        # 生成随机数矩阵 [n_simulations × 4]
        rand = np.random.random((n_simulations, len(p_locks)))
        
        # 判断每锁是否失效
        failures = rand < p_locks  # 布尔矩阵
        
        # 判断联合失效（四锁同时失效）
        joint_failures = np.all(failures, axis=1)
        
        # 统计
        n_joint = np.sum(joint_failures)
        
        if n_joint == 0:
            # 零失效：用95%置信上限
            estimate = 3.0 / n_simulations
        else:
            estimate = n_joint / n_simulations
        
        estimates.append(estimate)
    
    return np.array(estimates)

# 不同样本量的收敛分析
sample_sizes = np.logspace(3, 7, 15).astype(int)
results = []

for n in sample_sizes:
    estimates = monte_carlo_safety(p_locks, n, n_repeats=20)
    results.append({
        'n': n,
        'median': np.median(estimates),
        'lower': np.percentile(estimates, 5),
        'upper': np.percentile(estimates, 95)
    })

# ==================== 图1：单锁失效概率 ====================
fig1 = plt.figure(figsize=(8, 6))  # 设置尺寸为 8x6
ax1 = fig1.add_subplot(111)

colors = ['#FF6B6B', '#FFB347', '#6B9EFF', '#6BFF6B']
bars = ax1.bar(range(4), p_locks, color=colors, edgecolor='black', linewidth=1.5)

ax1.set_yscale('log')
ax1.set_ylim(1e-7, 5e-2)
ax1.set_xticks(range(4))
ax1.set_xticklabels(lock_names, fontsize=10)
ax1.set_ylabel('Failure Probability', fontsize=12)
ax1.set_title('Individual Safety Lock Failure Probability\n(Conservative Estimates)', fontsize=14)

# 添加数值标签
for i, (bar, p) in enumerate(zip(bars, p_locks)):
    height = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width()/2., height*1.5,
             f'$10^{{{int(np.log10(p))}}}$', 
             ha='center', va='bottom', fontsize=12, fontweight='bold')

# 调整文本框位置，避免与顶部数字重叠
textstr = 'Lock 1: Nutrient auxotrophy\nLock 2: Toxin-antitoxin\nLock 3: pH-sensitive lysis\nLock 4: Horizontal transfer block'
props = dict(boxstyle='round', facecolor='wheat', alpha=0.8)
ax1.text(0.02, 0.95, textstr, transform=ax1.transAxes, fontsize=10,
         verticalalignment='top', bbox=props)

plt.tight_layout()
plt.savefig('figure1_single_locks.png', dpi=600, bbox_inches='tight')
print("图1已保存为 figure1_single_locks.png (dpi=600)")

# ==================== 图2：联合失效概率收敛 ====================
fig2 = plt.figure(figsize=(8, 6))  # 与图1尺寸一致，均为 8x6
ax2 = fig2.add_subplot(111)

n_vals = [r['n'] for r in results]
medians = [r['median'] for r in results]
lowers = [r['lower'] for r in results]
uppers = [r['upper'] for r in results]

# 绘制置信区间
ax2.fill_between(n_vals, lowers, uppers, alpha=0.3, color='blue', label='95% CI')

# 绘制中位数
ax2.loglog(n_vals, medians, 'bo-', linewidth=2, markersize=6, label='Monte Carlo estimate')

# 理论线
ax2.axhline(y=p_theory, color='red', linestyle='--', linewidth=2, 
            label=f'Theoretical (independent): $10^{{{int(np.log10(p_theory))}}}$')

# FDA标准线
ax2.axhline(y=1e-12, color='green', linestyle=':', linewidth=2, alpha=0.7)
ax2.text(2e3, 2e-12, 'FDA Safety Standard\n($10^{-12}$)', fontsize=9, color='green')

# 零失效标注
ax2.annotate('No joint failure detected\nin $10^7$ simulations\n(Upper bound: $3\\times10^{-7}$)',
             xy=(1e6, 3e-7), xytext=(1e4, 1e-5),
             arrowprops=dict(arrowstyle='->', color='orange'),
             fontsize=10, color='orange',
             bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.8))

ax2.set_xlabel('Number of Simulations', fontsize=12)
ax2.set_ylabel('Estimated Joint Failure Probability', fontsize=12)
ax2.set_title('Monte Carlo Convergence Analysis\n(Four-Layer Safety System)', fontsize=14)
ax2.set_xlim(1e3, 2e7)
ax2.set_ylim(1e-20, 1e-2)
ax2.legend(loc='upper right', fontsize=10)
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('figure2_convergence.png', dpi=600, bbox_inches='tight')
print("图2已保存为 figure2_convergence.png (dpi=600)")

# ==================== 结果输出 ====================
print("\n" + "="*70)
print("Monte Carlo安全评估结果")
print("="*70)
print(f"最大模拟次数: {sample_sizes[-1]:,}")
print(f"重复次数: 20")
print(f"\n联合失效概率估计:")
print(f"  中位数: {medians[-1]:.2e}")
print(f"  95%置信区间: [{lowers[-1]:.2e}, {uppers[-1]:.2e}]")
print(f"\n理论值 (独立假设): {p_theory:.2e}")
print(f"与理论比值: {medians[-1]/p_theory:.2e}")
print(f"\n安全评估:")
print(f"  观测到联合失效次数: 0")
print(f"  95%置信上限: {3/sample_sizes[-1]:.2e}")
print(f"  FDA安全标准: 1e-12")
print(f"  结论: {'PASS' if uppers[-1] < 1e-12 else 'REVIEW'}")
print("="*70)