bw.py

def get_balloon_parameters(self):
    """get balloon model parameters."""
    alpha = 0.32
    E0 = 0.4
    eps = 0.5

    par = {
        "eps": eps,
        "E0": 0.4,
        "V0": 4.0,
        "alpha": alpha,
        "inv_alpha": 1.0 / alpha,
        "K1": 7.0 * E0,
        "K2": 2 * E0,
        "K3": 1 - eps,
        "taus": 1.54,
        "tauo": 0.98,
        "tauf": 1.44,
    }
    par["inv_tauo"] = 1.0 / par["tauo"]
    par["inv_taus"] = 1.0 / par["taus"]
    par["inv_tauf"] = 1.0 / par["tauf"]
    return par


def f_fmri(self, xin, x, t):

    E0 = self.E0
    xp = self.xp
    nn = self.num_nodes
    ns = self.num_sim
    inv_tauf = self.inv_tauf
    inv_tauo = self.inv_tauo
    inv_taus = self.inv_taus
    inv_alpha = self.inv_alpha

    dxdt = xp.zeros((4 * nn, ns)).astype(self.dtype)
    s = x[:nn, :]
    f = x[nn : 2 * nn, :]
    v = x[2 * nn : 3 * nn, :]
    q = x[3 * nn :, :]
    x0in = xin[:nn, :]  # use r

    dxdt[:nn, :] = x0in - inv_taus * s - inv_tauf * (f - 1.0)
    dxdt[nn : (2 * nn), :] = s
    dxdt[(2 * nn) : (3 * nn), :] = inv_tauo * (f - v ** (inv_alpha))
    dxdt[3 * nn :, :] = (inv_tauo) * (
        (f * (1.0 - (1.0 - E0) ** (1.0 / f)) / E0) - (v ** (inv_alpha)) * (q / v)
    )

    return dxdt


def integrate_fmri(self, yin, y, t):
    """
    Integrate Balloon model

    Parameters
    ----------
    yin : array [2*nn, ns]
        r and v time series, r is used as input
    y : array [4*nn, ns]
        state, update in place
    t : float
        time

    Returns
    -------
    yb : array [nn, ns]
        BOLD signal

    """

    V0 = self.V0
    K1 = self.K1
    K2 = self.K2
    K3 = self.K3

    nn = self.num_nodes
    self.heunDeterministic(yin, y, t)
    yb = V0 * (
        K1 * (1.0 - y[(3 * nn) :, :])
        + K2 * (1.0 - y[(3 * nn) :, :] / y[(2 * nn) : (3 * nn), :])
        + K3 * (1.0 - y[(2 * nn) : (3 * nn), :])
    )
    return yb