Source code for dynamo.vectorfield.Wang

from typing import Callable, Tuple

import numpy as np
from scipy import optimize
from scipy.optimize import OptimizeResult


[docs]def Wang_action(X_input: np.ndarray, F: Callable, D: float, dim: int, N: int, lamada_: float = 1) -> float: """Calculate action by path integral by Wang's method. Quantifying the Waddington landscape and biological paths for development and differentiation. Jin Wang, Kun Zhang, Li Xu, and Erkang Wang, PNAS, 2011 Args: X_input: The initial guess of the least action path. Default is a straight line connecting the starting and end path. F: The reconstructed vector field function. This is assumed to be time-independent. D: The diffusion constant. Note that this can be a space-dependent matrix. dim: The feature numbers of the input data. N: Number of waypoints along the least action path. lamada_: Regularization parameter Returns: The action function calculated by the Hamilton-Jacobian method. """ X_input = X_input.reshape((int(dim), -1)) if len(X_input.shape) == 1 else X_input delta, delta_l = delta_delta_l(X_input) V_m = np.zeros((N, 1)) F_l = np.zeros((N, 1)) E_eff = np.zeros((N, 1)) for i in range(N - 1): F_m = F(X_input[:, i]).reshape((1, -1)) V_m[i] = V(F, D, X_input[:, i]) E_eff[i] = np.sum(F(X_input[:, i]) ** 2) / (4 * D) - V_m[i] F_l[i] = F_m.dot(delta[:, i]) / delta_l[i] P = np.sum((delta_l - np.linalg.norm(X_input[:, N - 1] - X_input[:, 0]) / N) ** 2) S_HJ = np.sum((np.sqrt((E_eff + V_m[:N]) / D) - 1 / (2 * D) * F_l) * delta_l) + lamada_ * P print(S_HJ) return S_HJ
def V_jacobina(F, X): import numdifftools as nda V_jacobina = nda.Jacobian(F) return V_jacobina(X) def V(F: Callable, D: float, X: np.ndarray) -> np.ndarray: """Calculate V Args: F: `Function` The reconstructed vector field function D: `float` The diffusion constant X: `nummpy.ndarray` The input coordinates corresponding to the cell states. Returns: Returns V """ V = 1 / (4 * D) * np.sum(F(X) ** 2) + 1 / 2 * np.trace(V_jacobina(F, X)) return V def delta_delta_l(X_input) -> Tuple[np.ndarray, float]: """Calculate delta_L""" delta = np.diff(X_input, 1, 1) delta_l = np.sqrt(np.sum(delta**2, 0)) return delta, delta_l
[docs]def Wang_LAP( F: Callable, n_points: int, point_start: np.ndarray, point_end: np.ndarray, D: float = 0.1, lambda_: float = 1 ) -> OptimizeResult: """Calculating least action path based methods from Jin Wang and colleagues (http://www.pnas.org/cgi/doi/10.1073/pnas.1017017108) Args: F: The reconstructed vector field function n_points: The number of points along the least action path. point_start: The matrix for storing the coordinates (gene expression configuration) of the start point (initial cell state). point_end: The matrix for storing the coordinates (gene expression configuration) of the end point (terminal cell state). D: The diffusion constant. Note that this can be a space-dependent matrix. lamada_: Regularization parameter Returns: The least action path and the action way of the inferred path. """ initpath = point_start.dot(np.ones((1, n_points + 1))) + (point_end - point_start).dot( np.linspace(0, 1, n_points + 1, endpoint=True).reshape(1, -1) ) dim, N = initpath.shape # update this optimization method res = optimize.basinhopping( Wang_action, x0=initpath, minimizer_kwargs={"args": (F, D, dim, N, lambda_)}, ) return res
[docs]def transition_rate(X_input: np.ndarray, F: Callable, D: float = 0.1, lambda_: float = 1) -> np.ndarray: """Calculate the rate to convert from one cell state to another cell state by taking the optimal path. In the small noise limit (D -> 0) the Wentzell-Freidlin theory states that the transition rate from one basin to another one to a leading order is related to the minimum action corresponding zero energy path (Eeff = 0) connecting the starting fixed point and the saddle point x_{sd} by k \approx exp(−S0 (x_{sd})). To take into account that for finite noise, the actual optimal path bypasses the saddle point, in Eqn. 2 of the main text a transition rate is actually estimated by the action of the whole path connecting the two fixed points, giving that the portion of the path following the vector field contributes zero action. Here we have neglected some pre-exponential factor (see Eq. 5.24 of reference [15]), which is expected to be on the order of 1 [12]. (Reference: Epigenetic state network approach for describing cell phenotypic transitions. Ping Wang, Chaoming Song, Hang Zhang, Zhanghan Wu, Xiao-Jun Tian and Jianhua Xing) Args: X_input: The initial guess of the least action path. Default is a straight line connecting the starting and end path. F: The reconstructed vector field function D: The diffusion constant. Note that this can be a space-dependent matrix. lamada_: Regularization parameter Returns: The transition to convert from one cell state to another. """ res = Wang_LAP(X_input, F, D=D, lambda_=lambda_) r = np.exp(-res) return r
[docs]def MFPT(X_input: np.ndarray, F: Callable, D: float = 0.1, lambda_: float = 1) -> float: """Calculate the MFPT (mean first passage time) to convert from one cell state to another cell state by taking the optimal path. The mean first-passage time (MFPT) defines an average timescale for a stochastic event to first occur. The MFPT maps a multi-step kinetic process to a coarse-grained timescale for reaching a final state, having started at some initial state. The inverse of the MFPT is an effective rate of the overall reaction. (reference: Mean First-Passage Times in Biology Nicholas F. Polizzi,a Michael J. Therien,b and David N. Beratan) Args: X_input: The initial guess of the least action path. Default is a straight line connecting the starting and end path. F: The reconstructed vector field function D: The diffusion constant. Note that this can be a space-dependent matrix. lamada_: Regularization parameter Returns: The transition to convert from one cell state to another. """ r = transition_rate(X_input, F, D=D, lambda_=lambda_) t = 1 / r return t