eqtk.solve¶
- eqtk.solve(c0, N=None, K=None, logK=None, A=None, G=None, names=None, units=None, G_units=None, solvent_density=None, T=293.15, return_log=False, normal_A=True, tol=1e-07, tol_zero=1e-12, max_iters=1000, delta_bar=1000.0, eta=0.125, min_delta=1e-12, max_trials=100, perturb_scale=100.0)¶
Solve for equilibrium concentrations of all species in a dilute solution.
- Parameters
c0 (array_like, dict, Series, or DataFrame, shape (n_points, n_compounds) or (n_compounds, )) – Each row contains the total “initial” concentration of all possible chemical species in solution. The equilibrium concentration of all species is computed for each row in c0. c0[i, j] is the initial concentration of compound j for calculation i. c0 may also be passed as a Pandas Series where the indices contain the name of the chemical species and each value is the “initial concentration.” c0 may also be passed as a Pandas DataFrame where each row contains the total “initial” concentration of all possible compounds in solution and the column names are the names of the chemical species. If c0 is passed as a dict, the dict must be convertible to a Pandas Series or DataFrame as pd.Series(c0) or pd.DataFrame(c0).
N (array_like or DataFrame, default None) – Stoichiometic matrix. N[r, j] = the stoichiometric coefficient of compound j in chemical reaction r. All rows of N must be linearly independent. If entered as a DataFrame, the name of chemical species j is N.columns[j]. Optionally, column ‘equilibrium constant’ contains the equilibrium constants for each reaction in units commensurate with those of c0. If N is given, A and G cannot be given.
K (array_like, shape (n_reactions,), default None) – K[r] is the equilibrium constant for chemical reaction r in units commensurate with those of c0. If N is given as a DataFrame with an ‘equilibrium constant’ column, K should not be supplied. If K is given, A and G cannot be given.
logK (array_like, shape (n_reactions,), default None) – logK[r] is the natural logarithm of the equilibrium constant for chemical reaction r. If logK is specified, the concentrations must all be dimensionless (units=None). If N is given as a DataFrame with a ‘log equilibrium constant’ column, logK should not be supplied. If K is given, A, G, and K cannot be given.
A (array_like or DataFrame, n_compounds columns) – Conservation matrix. If c is the output, then A @ c0 = A @ c. All entries must be nonnegative and the rows of A must be linearly independent. If entered as a DataFrame, the name of chemical species j is A.columns[j]. If A is given, G must be given, and N and K cannot be given.
G (array_like, shape (n_compounds, ), default None) – G[j] is the free energy of chemical species j in units specified by G_units. If G is given, A must be given, and N and K cannot be given.
names (list or tuple of str, default None, optional) – The names of the chemical species. Names are inferred if N or A is given as a DataFrame, in which case names is unnecessary.
units (string or None, default None) – The units of the concentrations inputted as c0. The output is also in these units. Allowable values are {None, ‘mole fraction’, ‘molar’, ‘M’, ‘millimolar’, ‘mM’, ‘micromolar’, ‘uM’, ‘µM’, ‘nanomolar’, ‘nM’, ‘picomolar’, ‘pM’}. If None, concentrations are considered to be dimensionless. The equilibrium constants given by K must have corresponding units.
G_units (string, default None) – Units in which free energy is given. If None or ‘kT’, the free energies are specified in units of of the thermal energy kT. Allowable values are {None, ‘kT’, kcal/mol’, ‘J’, ‘J/mol’, ‘kJ/mol’, ‘pN-nm’}.
solvent_density (float, default None) – The density of the solvent in units commensurate with the units keyword argument. Default (None) assumes the solvent is water, and its density is computed at the temperature specified by the T keyword argument.
T (float, default = 293.15) – Temperature, in Kelvin, of the solution. When N and K are given, T is ignored if solvent_density is given or if units is None. If A and G are given, T is ignored when units and G_units are both None.
return_log (bool, default False) – If True, return the natural logarithm of the concentrations.
normal_A (bool, default True) – If True, perform a transformation on A such that its rows are orthonormal. If False, use inputted A directly as the conservation matrix. This is ignored if A is not specified; the resulting conservation matrix in that case has orthonormal rows by construction.
tol (float, default 1e-7) – Tolerance for convergence. The absolute tolerance for a given initial concentration c0 (a one-dimensional array) for the conservation conditions are tol * np.dot(A, c0) except when an entry in np.dot(A, c0) is zero. If that is the case for entry i, then the absolute tolerance is tol * np.max(A[i] * c0). If all entries in A[i] * c0 are zero, which only happens with c0 = 0, the absolute tolerance is tol_zero.
tol_zero (float, default 1e-12) – Absolute tolerance for convergence when the initial concentrations are all zero. This cannot really be set a priori; it depends on the scale of the dimensionless concentrations. By default, assume an absolute tolerance consistent with tol and millimolar mole fractions.
max_iters (int, default 1000) – Maximum number of iterations allowed in trust region method.
delta_bar (float, default 1000.0) – Maximum step size allowed in the trust region method.
eta (float, default 0.125) – Value for eta in the trust region method. eta must satisfy 0 < eta < 0.25.
min_delta (float, default 1e-12) – Minimal allowed radius of the trust region. When the trust region radius gets below min_delta, the trust region iterations stop, and a final set of Newton steps is attempted.
max_trials (int, default 100) – In the event that an attempt to solve does not converge, the solver tries again with different initial guesses. This continues until max_trials failures.
perturb_scale (float, default 100.0) – Multiplier on random perturbations to the initial guesses as new ones are generated.
- Returns
c – Equilibrium concentrations of all species. c[i, j] is the equilibrium concentration of species j for initial concentrations given by c0[i, :] in units given by units. If c0 is inputted as a DataFrame or names is not None, then c is a DataFrame with columns given by names or with the same columns (without ‘equilibrium constant’) as c0. Otherwise, c is returned as a Numpy array with the same shape as c0. If return_log is True, then the return value is the natural logarithm of the dimensionless concentrations.
- Return type
array or DataFrame, shape c0.shape
- Raises
ValueError – If input is in any way invalid
RuntimeError – If the trust region algorithm failed to converge
Notes
Uses an elliptical trust region optimization to find the equilibrium concentrations. See 1 for algorithmic details, as well as definitions of the parameters associated with the trust region algorithm.
In practice, the trust region parameters should not be adjusted from their default values.
References
- 1
Nocedal and Wright, Numerical Optimization, Second Edition, Springer, 2006, Chapter 4.
Examples
1) Find the equilibrium concentrations of a solution containing species A, B, C, AB, BB, and BC that can undergo chemical reactions
A ⇌ C, K = 50 (dimensionless)
A + C ⇌ AB K = 10 (1/mM)
B + B ⇌ BB K = 40 (1/mM)
B + C ⇌ BC K = 100 (1/mM)
with initial concentrations of [A]₀ = 1.0 mM, [B]₀ = 3.0 mM.
>>> N = np.array([[-1, 0, 1, 0, 0, 0], ... [ 1, 1, 0, -1, 0, 0], ... [ 0, 2, 0, 0, -1, 0], ... [ 0, 1, 1, 0, 0, -1]]) >>> K = np.array([50.0, 0.1, 0.025, 0.01]) >>> c0 = np.array([1.0, 3.0, 0.0, 0.0, 0.0, 0.0]) >>> eqtk.solve(c0, N=N, K=K, units='mM') array([0.00121271, 0.15441164, 0.06063529, 0.00187256, 0.95371818, 0.93627945])
2) Compute a titration curve for the same reaction system as in example (1) with [A]₀ = 1.0 mM. Consider B being titrated from [B]₀ = 0.0 to 3.0 and only use four titration points for display purposes.
>>> names = ['A', 'B', 'C', 'AB', 'BB', 'BC'] >>> df_N = pd.DataFrame( ... data=[[-1, 0, 1, 0, 0, 0], ... [ 1, 1, 0, -1, 0, 0], ... [ 0, 2, 0, 0, -1, 0], ... [ 0, 1, 1, 0, 0, -1]], ... columns=names ... ) >>> df_N['equilibrium constant'] = [50.0, 0.1, 0.025, 0.01] >>> df_c0 = pd.DataFrame(data=np.zeros((4, 6)), columns=names) >>> df_c0['A'] = 1.0 >>> df_c0['B'] = np.arange(4) >>> df_c = eqtk.solve(df_c0, N=df_N, units='mM') >>> df_c.loc[:, ~df_c.columns.str.contains('__0')] # Don't disp. c0 A B C AB BB BC 0 0.019608 0.000000 0.980392 0.000000 0.000000 0.000000 1 0.003750 0.043044 0.187513 0.001614 0.074110 0.807122 2 0.001656 0.110347 0.082804 0.001827 0.487057 0.913713 3 0.001213 0.154412 0.060635 0.001873 0.953718 0.936279
3) Find the equilibrium concentrations of a solution containing species A, B, C, AB, and AC that have free energies (in units of kT).
A : 0.0
B : 1.0
C : -2.0
AB : -3.0
AC : -7.0,
with initial mole fractions x0_A = 0.005, x0_B = 0.001, and x0_C = 0.002. The ordering of the compounds in the example is A, B, C, AB, AC.
>>> A = np.array([[ 1, 0, 0, 1, 1], ... [ 0, 1, 0, 1, 0], ... [ 0, 0, 1, 0, 1]]) >>> G = np.array([0.0, 1.0, -2.0, -3.0, -7.0]) >>> c0 = np.array([0.005, 0.001, 0.002, 0.0, 0.0]) >>> eqtk.solve(c0, A=A, G=G) array([0.00406569, 0.00081834, 0.00124735, 0.00018166, 0.00075265])
4) Competition assay. Species A binds both B and C with specified dissociation constants, Kd.
AB ⇌ A + B, Kd = 50 nM
AC ⇌ A + C Kd = 10 nM
Initial concentrations of [A]_0 = 2.0 uM, [B]_0 = 0.05 uM, [C]_0 = 1.0 uM. The ordering of the compounds in the example is A, B, C, AB, AC.
>>> rxns = ''' ... AB <=> A + B ; 0.05 ... AC <=> A + C ; 0.01''' >>> N = eqtk.parse_rxns(rxns) >>> c0 = pd.Series({'A': 2.0, 'B': 0.05, 'C': 1, 'AB': 0, 'AC': 0}) >>> eqtk.solve(c0, N=N, units='µM') A__0 2.000000 B__0 0.050000 C__0 1.000000 AB__0 0.000000 AC__0 0.000000 A 0.962749 B 0.002469 C 0.010280 AB 0.047531 AC 0.989720 dtype: float64