def METHOD_NAME(self, rates, max_time):
    """ Returns raw times, mol counts when reactions occur. """
    if len(rates) != self.n_parameters():
        raise ValueError(
            'This model should have ' + str(self.n_parameters())
            + ' parameter(s).')
    # Setting the current propensities and summing them up
    current_propensities = self._propensities(self._x0, rates)
    prop_sum = np.sum(current_propensities)
    # Initial time and count
    t = 0
    x = np.array(self._x0)
    # Run Gillespie SSA, calculating time until next reaction, deciding
    # which reaction, and applying it
    mol_count = [np.array(x)]
    time = [t]
    while prop_sum > 0 and t <= max_time:
        r_1, r_2 = np.random.uniform(0, 1), np.random.uniform(0, 1)
        t += -np.log(r_1) / prop_sum
        s = 0
        r = 0
        while s <= r_2 * prop_sum:
            s += current_propensities[r]
            r += 1
        x += self._V[r - 1]
        # Calculate new current propensities
        current_propensities = self._propensities(x, rates)
        prop_sum = np.sum(current_propensities)
        # Store new values
        time.append(t)
        mol_count.append(np.copy(x))
    return np.array(time), np.array(mol_count)