NeuronSynthesis/model.py at main · fcavarretta-ualr/NeuronSynthesis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import numpy as np

def mean_Z(x, E_z0, beta, alpha):
    """
    Compute the mean of Z(x) using kappa = beta - alpha.

    Parameters:
        x (float or array): Position(s)
        E_z0 (float): Expected value of z0
        beta (float): Growth rate parameter
        alpha (float): Decay rate parameter

    Returns:
        float or np.ndarray: E[Z(x)]
    """
    kappa = beta - alpha
    return E_z0 * np.exp(kappa * x)

def var_Z(x, E_z0, Var_z0, beta, alpha):
    """
    Compute the variance of Z(x) using kappa = beta - alpha.

    Parameters:
        x (float or array): Position(s)
        E_z0 (float): Expected value of z0
        Var_z0 (float): Variance of z0
        beta (float): Growth rate parameter
        alpha (float): Decay rate parameter

    Returns:
        float or np.ndarray: Var[Z(x)]
    """
    kappa = beta - alpha

    if kappa != 0:
        exp_kx = np.exp(kappa * x)
        term1 = E_z0 * (2 * beta - kappa) * exp_kx * (exp_kx - 1) / kappa
        term2 = Var_z0 * np.exp(2 * kappa * x)
        return term1 + term2
    else:
        return 2 * E_z0 * beta * x + Var_z0


def cov_Z(x, E_z0, Var_z0, beta, alpha):
    """
    Compute the variance of Z(x) using kappa = beta - alpha.

    Parameters:
        x (float or array): Position(s)
        E_z0 (float): Expected value of z0
        Var_z0 (float): Variance of z0
        beta (float): Growth rate parameter
        alpha (float): Decay rate parameter

    Returns:
        float or np.ndarray: Var[Z(x)]
    """
    kappa = beta - alpha

    return Var_z0 * np.exp(kappa * x)


def mean_B(x, E_z0, beta, alpha):
    """
    Compute the mean number of bifurcations B(x).

    Parameters:
        x (float or array): Position(s)
        E_z0 (float): Expected initial number of branches z0
        beta (float): Growth rate
        alpha (float): Pruning or death rate

    Returns:
        float or np.ndarray: E[B(x)]
    """
    kappa = beta - alpha
    if abs(kappa) > 1e-5:
        return beta * E_z0 * (np.exp(kappa * x) - 1) / kappa
    else:
        return beta * E_z0 * x


def var_B(x, E_z0, Var_z0, beta, alpha):
    """
    Compute the variance of the number of bifurcations B(x).

    Parameters:
        x (float or array): Position(s)
        E_z0 (float): Expected initial number of branches z0
        Var_z0 (float): Variance of z0
        beta (float): Growth rate
        alpha (float): Pruning or death rate

    Returns:
        float or np.ndarray: Var[B(x)]
    """

    kappa = beta - alpha

    if abs(kappa) > 1e-5:
        exp_kappa = np.exp(kappa * x)
        exp_2kappa = np.exp(2 * kappa * x)

        term1 = beta * E_z0 * (exp_kappa - 1) / kappa
        term2 = beta**2 * E_z0 * (2 * beta - kappa) * (exp_2kappa - 2 * kappa * x * exp_kappa - 1) / kappa**3
        term3 = beta**2 * Var_z0 * ((exp_kappa - 1) / kappa)**2

        return term1 + term2 + term3
    else:
        return ( beta * x ) ** 2 * (2 / 3 * beta * E_z0 * x + Var_z0) + beta * E_z0 * x