function dxdt = sys(t,X,Q)
    [be, delta, d1, k2, lambda1, Kb] = Q
    d2 = 0.01
    k1 = 8e-7
    lambda2 = 31.98
    eps = .1
    de = 0.25
    Nt = 100
    f = 0.34
    m1 = m2 = 1e-5
    lambda_e = 1
    delta_e = 0.1
    p1 = p2 = 1
    Kd = 1
    T1 = X(1)
    T2 = X(2)
    T1s = X(3)
    T2s = X(4)
    V = X(5)
    E = X(6)

    dT1 = lambda1 - d1*T1 - (1-eps)*k1*V*T1
    dT2 = lambda2 - d2*T2 - (1-f*eps)*k2*V*T2
    dT1s = (1-eps)*k1*V*T1 - delta*T1s - m1*E*T1s
    dT2s = (1 - f*eps)*k2*V*T2 - delta*T2s - m2*E*T2s
    dV = Nt*delta*(T1s + T2s) - c*V -((1-eps)*p1*k1*T1 + (1-f*eps)*p2*k2*T2)*V
    dE = lambda_e + E*(be*(T1s + T2s))/(T1s + T2s + Kb) - E*(de*(T1s + T2s))/(T1s + T2s + Kd) - delta_e*E

    dxdt = [dT1,dT2,dT1s,dT2s,dV,dE]