def reactions(mu, states):
    """Executes Michaelis-Menten chemical reactions.

    :Arguments:
        mu : int
            Index of the equation
        states : NumPy array
            Current states of the system (E, S, ES, P)

    :Returns:
        NumPy array : Updated state vector (E, S, ES, P)
    """
    if mu == 1:
        # E + S -> ES
        if states[0] > 0 and states[1] > 0:
            states[0] -= 1
            states[1] -= 1
            states[2] += 1
    elif mu == 2:
        # ES -> E + S
        if states[2] > 0:
            states[0] += 1
            states[1] += 1
            states[2] -= 1
    elif mu == 3:
        # ES -> E + P
        if states[2] > 0:
            states[0] += 1
            states[2] -= 1
            states[3] += 1
    else:
        raise ValueError('Reaction mu = %d does not exist.' % mu)

    return states