import logging; logger = logging.getLogger('matplotlib'); logger.setLevel(logging.INFO) # or notebook filled with logging from collections import OrderedDict, defaultdict import math import re import time from IPython.display import Image, Latex, display import matplotlib.pyplot as plt import sympy from pyodesys.symbolic import ScaledSys from pyodesys.native.cvode import NativeCvodeSys from chempy import Substance, Equilibrium, Reaction, ReactionSystem from chempy.kinetics.ode import get_odesys from chempy.kinetics.rates import MassAction from chempy.printing.tables import UnimolecularTable, BimolecularTable from chempy.thermodynamics.expressions import EqExpr from chempy.util.graph import rsys2graph from chempy.util.pyutil import defaultkeydict %matplotlib inline substances = OrderedDict([ ('N', Substance('N', composition={'protein': 1}, latex_name='[N]')), ('U', Substance('U', composition={'protein': 1}, latex_name='[U]')), ('A', Substance('A', composition={'protein': 1}, latex_name='[A]')), ('L', Substance('L', composition={'ligand': 1}, latex_name='[L]')), ('NL', Substance('NL', composition={'protein': 1, 'ligand': 1}, latex_name='[NL]')), ]) def _gibbs(args, T, R, backend, **kwargs): H, S, Cp, Tref = args H2 = H + Cp*(T - Tref) S2 = S + Cp*backend.log(T/Tref) return backend.exp(-(H2 - T*S2)/(R*T)) def _eyring(args, T, R, k_B, h, backend, **kwargs): H, S = args return k_B/h*T*backend.exp(-(H - T*S)/(R*T)) Gibbs = EqExpr.from_callback(_gibbs, parameter_keys=('temperature', 'R'), argument_names=('H', 'S', 'Cp', 'Tref')) Eyring = MassAction.from_callback(_eyring, parameter_keys=('temperature', 'R', 'k_B', 'h'), argument_names=('H', 'S')) thermo_dis = Gibbs(unique_keys=('He_dis', 'Se_dis', 'Cp_dis', 'Tref_dis')) thermo_u = Gibbs(unique_keys=('He_u', 'Se_u', 'Cp_u', 'Tref_u')) # ([He_u_R, Se_u_R, Cp_u_R, Tref]) kinetics_agg = Eyring(unique_keys=('Ha_agg', 'Sa_agg')) # EyringMassAction([Ha_agg, Sa_agg]) kinetics_as = Eyring(unique_keys=('Ha_as', 'Sa_as')) kinetics_f = Eyring(unique_keys=('Ha_f', 'Sa_f')) eq_dis = Equilibrium({'NL'}, {'N', 'L'}, thermo_dis, name='ligand-protein dissociation') eq_u = Equilibrium({'N'}, {'U'}, thermo_u, {'L'}, {'L'}, name='protein unfolding') r_agg = Reaction({'U'}, {'A'}, kinetics_agg, {'L'}, {'L'}, name='protein aggregation') rsys = ReactionSystem( eq_dis.as_reactions(kb=kinetics_as, new_name='ligand-protein association') + eq_u.as_reactions(kb=kinetics_f, new_name='protein folding') + (r_agg,), substances, name='4-state CETSA system') vecs, comp = rsys.composition_balance_vectors() names = rsys.substance_names() dict(zip(comp, [dict(zip(names, v)) for v in vecs])) rsys2graph(rsys, '4state.png', save='.', include_inactive=False) Image('4state.png') rsys uni, not_uni = UnimolecularTable.from_ReactionSystem(rsys) bi, not_bi = BimolecularTable.from_ReactionSystem(rsys) assert not (not_bi & not_uni), "Only uni- & bi-molecular reactions expected" uni bi def pretty_replace(s, subs=None): if subs is None: subs = { 'Ha_(\w+)': r'\\Delta_{\1}H^{\\neq}', 'Sa_(\w+)': r'\\Delta_{\1}S^{\\neq}', 'He_(\w+)': r'\\Delta_{\1}H^\\circ', 'Se_(\w+)': r'\\Delta_{\1}S^\\circ', 'Cp_(\w+)': r'\\Delta_{\1}\,C_p', 'Tref_(\w+)': r'T^{\\circ}_{\1}', } for pattern, repl in subs.items(): s = re.sub(pattern, repl, s) return s def mk_Symbol(key): if key in substances: arg = substances[key].latex_name else: arg = pretty_replace(key.replace('temperature', 'T')) return sympy.Symbol(arg) autosymbols = defaultkeydict(mk_Symbol) rnames = {} for rxn in rsys.rxns: rnames[rxn.name] = rxn.name.replace(' ', '~').replace('-','-') rate_expr_str = sympy.latex(rxn.rate_expr()(autosymbols, backend=sympy, reaction=rxn)) lstr = r'$r(\mathrm{%s}) = %s$' % (rnames[rxn.name], rate_expr_str) display(Latex(lstr)) ratexs = [autosymbols['r(\mathrm{%s})' % rnames[rxn.name]] for rxn in rsys.rxns] rates = rsys.rates(autosymbols, backend=sympy, ratexs=ratexs) for k, v in rates.items(): display(Latex(r'$\frac{[%s]}{dt} = %s$' % (k, sympy.latex(v)))) default_c0 = defaultdict(float, {'N': 1e-9, 'L': 1e-8}) params = dict( R=8.314472, # or N_A & k_B k_B=1.3806504e-23, h=6.62606896e-34, # k_B/h == 2.083664399411865e10 K**-1 * s**-1 He_dis=-45e3, Se_dis=-400, Cp_dis=1.78e3, Tref_dis=298.15, He_u=60e3, Cp_u=20.5e3, Tref_u=298.15, Ha_agg=106e3, Sa_agg=70, Ha_as=4e3, Sa_as=-10, Ha_f=90e3, Sa_f=50, temperature=50 + 273.15 ) def Se0_from_Tm(Tm, token): dH0, T0, dCp = params['He_'+token], params['Tref_'+token], params['Cp_'+token] return dH0/Tm + (Tm-T0)*dCp/Tm - dCp*math.log(Tm/T0) params['Se_u'] = Se0_from_Tm(48.2+273.15, 'u') params['Se_u'] params_c0 = default_c0.copy() params_c0.update(params) for rxn in rsys.rxns: print('%s: %.5g' % (rxn.name, rxn.rate_expr()(params_c0, reaction=rxn))) odesys, extra = get_odesys(rsys, include_params=False, SymbolicSys=ScaledSys, dep_scaling=1e9) len(odesys.exprs) # how many (symbolic) expressions are there in this representation? h0max = extra['max_euler_step_cb'](0, default_c0, params) h0max def integrate_and_plot(system, c0=None, first_step=None, t0=0, stiffness=False, nsteps=9000, **kwargs): if c0 is None: c0 = default_c0 if first_step is None: first_step = h0max*1e-11 tend = 3600*24 t_py = time.time() kwargs['atol'] = kwargs.get('atol', 1e-11) kwargs['rtol'] = kwargs.get('rtol', 1e-11) res = system.integrate([t0, tend], c0, params, integrator='cvode', nsteps=nsteps, first_step=first_step, **kwargs) t_py = time.time() - t_py if stiffness: plt.subplot(1, 2, 1) _ = system.plot_result(xscale='log', yscale='log') _ = plt.legend(loc='best') plt.gca().set_ylim([1e-16, 1e-7]) plt.gca().set_xlim([1e-11, tend]) if stiffness: if stiffness is True: stiffness = 0 ratios = odesys.stiffness() plt.subplot(1, 2, 2) plt.yscale('linear') plt.plot(odesys._internal[0][stiffness:], ratios[stiffness:]) for k in ('time_wall', 'time_cpu'): print('%s = %.3g' % (k, res[2][k]), end=', ') print('time_python = %.3g' % t_py) return res _, _, info = integrate_and_plot(odesys) assert info['internal_yout'].shape[1] == 5 {k: v for k, v in info.items() if not k.startswith('internal')} native = NativeCvodeSys.from_other(odesys, first_step_expr=0*odesys.indep) _, _, info_native = integrate_and_plot(native) {k: v for k, v in info_native.items() if not k.startswith('internal')} info['time_wall']/info_native['time_wall'] from chempy.kinetics._native import get_native native2 = get_native(rsys, odesys, 'cvode') _, _, info_native2 = integrate_and_plot(native2, first_step=0.0) {k: v for k, v in info_native2.items() if not k.startswith('internal')} cses, (jac_in_cse,) = odesys.be.cse(odesys.get_jac()) jac_in_cse odesys.jacobian_singular() A, comp_names = rsys.composition_balance_vectors() A, comp_names, list(rsys.substances.keys()) y0 = {odesys[k]: sympy.Symbol(k+'0') for k in rsys.substances.keys()} analytic_L_N = extra['linear_dependencies'](['L', 'N']) analytic_L_N(None, y0, None, sympy) assert len(analytic_L_N(None, y0, None, sympy)) > 0 # ensure the callback is idempotent analytic_L_N(None, y0, None, sympy), list(enumerate(odesys.names)) from pyodesys.symbolic import PartiallySolvedSystem no_invar = dict(linear_invariants=None, linear_invariant_names=None) psysLN = PartiallySolvedSystem(odesys, analytic_L_N, **no_invar) print(psysLN.be.cse(psysLN.get_jac())[1][0]) psysLN['L'], psysLN.jacobian_singular(), len(psysLN.exprs) psysLA = PartiallySolvedSystem(odesys, extra['linear_dependencies'](['L', 'A']), **no_invar) print(psysLA.be.cse(psysLA.get_jac())[1][0]) psysLA['L'], psysLA.jacobian_singular() plt.figure(figsize=(12,4)) plt.subplot(1, 2, 1) _, _, info_LN = integrate_and_plot(psysLN, first_step=0.0) assert info_LN['internal_yout'].shape[1] == 3 plt.subplot(1, 2, 2) _, _, info_LA = integrate_and_plot(psysLA, first_step=0.0) assert info_LA['internal_yout'].shape[1] == 3 ({k: v for k, v in info_LN.items() if not k.startswith('internal')}, {k: v for k, v in info_LA.items() if not k.startswith('internal')}) from pyodesys.symbolic import SymbolicSys psys_root = SymbolicSys.from_other(psysLN, roots=[psysLN['N'] - psysLN['A']]) psys_root.roots psysLN['N'] psysLN.analytic_exprs psysLN.names psysLN.dep tout1, Cout1, info_root = integrate_and_plot(psys_root, first_step=0.0, return_on_root=True) print('Time at which concnetrations of N & A are equal: %.4g' % (tout1[-1])) xout2, yout2, info_LA = integrate_and_plot(psysLA, first_step=0.0, t0=tout1[-1], c0=dict(zip(odesys.names, Cout1[-1, :]))) print('\troot\tLA\troot+LA\tLN') for k in 'n_steps nfev njev'.split(): print('\t'.join(map(str, (k, info_root[k], info_LA[k], info_root[k] + info_LA[k], info_LN[k])))) from pyodesys.symbolic import symmetricsys logexp = lambda x: sympy.log(x + 1e-20), lambda x: sympy.exp(x) - 1e-20 def psimp(exprs): return [sympy.powsimp(expr.expand(), force=True) for expr in exprs] LogLogSys = symmetricsys(logexp, logexp, exprs_process_cb=psimp) unscaled_odesys, unscaled_extra = get_odesys(rsys, include_params=False) tsys = LogLogSys.from_other(unscaled_odesys) unscaledLN = PartiallySolvedSystem(unscaled_odesys, unscaled_extra['linear_dependencies'](['L', 'N']), **no_invar) unscaledLA = PartiallySolvedSystem(unscaled_odesys, unscaled_extra['linear_dependencies'](['L', 'A']), **no_invar) assert sorted(unscaledLN.free_names) == sorted(['U', 'A', 'NL']) assert sorted(unscaledLA.free_names) == sorted(['U', 'N', 'NL']) tsysLN = LogLogSys.from_other(unscaledLN) tsysLA = LogLogSys.from_other(unscaledLA) _, _, info_t = integrate_and_plot(tsys, first_step=0.0) {k: info_t[k] for k in ('nfev', 'njev', 'n_steps')} native_tLN = NativeCvodeSys.from_other(tsysLN) _, _, info_tLN = integrate_and_plot(native_tLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-9) {k: info_tLN[k] for k in ('nfev', 'njev', 'n_steps')} _, _, info_tLN = integrate_and_plot(tsysLN, first_step=1e-9, nsteps=18000, atol=1e-8, rtol=1e-8) {k: info_tLN[k] for k in ('nfev', 'njev', 'n_steps')} _, _, info_tLA = integrate_and_plot(tsysLA, first_step=0.0) {k: info_tLA[k] for k in ('nfev', 'njev', 'n_steps')} print(open(next(filter(lambda s: s.endswith('.cpp'), native2._native._written_files))).read())