How to use the chempy.kinetics.rates.MassAction function in chempy

def _reg_unique(expr, rxn=None):
        if not isinstance(expr, Expr):
            raise NotImplementedError("Currently only Expr sub classes are supported.")

        if isinstance(expr, MassAction):
            if expr.args is None:
                uk, = expr.unique_keys
                if uk not in substitutions:
                    unique[uk] = None
                    _reg_unique_unit(uk, _get_arg_dim(expr, rxn), 0)
                    return
            else:
                arg, = expr.args
                if isinstance(arg, Symbol):
                    uk, = arg.unique_keys
                    if uk not in substitutions:
                        unique[uk] = None
                        _reg_unique_unit(uk, _get_arg_dim(expr, rxn), 0)
                        return

        if expr.args is None:

def get_rsys(kf, kb):
    rf = MassAction([kf], unique_keys=['kf'])
    rb = MassAction([kb], unique_keys=['kb'])
    fw = Reaction({'Fe+3', 'SCN-'}, {'FeSCN+2'}, rf)
    bw = Reaction({'FeSCN+2'}, {'Fe+3', 'SCN-'}, rb)
    return ReactionSystem([fw, bw], 'Fe+3 SCN- FeSCN+2'.split())

def get_rsys(kf, kb):
    rf = MassAction([kf], unique_keys=['kf'])
    rb = MassAction([kb], unique_keys=['kb'])
    fw = Reaction({'Fe+3', 'SCN-'}, {'FeSCN+2'}, rf)
    bw = Reaction({'FeSCN+2'}, {'Fe+3', 'SCN-'}, rb)
    return ReactionSystem([fw, bw], 'Fe+3 SCN- FeSCN+2'.split())

if globals_ is None:
        globals_ = get_parsing_context()
    parts = line.rstrip('\n').split(';')
    stoich = parts[0].strip()
    if len(parts) > 2:
        kwargs.update(eval('dict('+';'.join(parts[2:])+'\n)', globals_ or {}))
    if len(parts) > 1:
        param = parts[1].strip()
    else:
        param = kwargs.pop('param', 'None')

    if isinstance(param, str):
        if param.startswith("'") and param.endswith("'") and "'" not in param[1:-1]:
            from ..kinetics.rates import MassAction
            from ._expr import Symbol
            param = MassAction(Symbol(unique_keys=(param[1:-1],)))
        else:
            param = None if globals_ is False else eval(param, globals_)

    if token not in stoich:
        raise ValueError("Missing token: %s" % token)

    reac_prod = [[y.strip() for y in x.split(' + ')] for x in stoich.split(token)]

    act, inact = [], []
    for elements in reac_prod:
        act.append(_parse_multiplicity([x for x in elements if not x.startswith('(')], substance_keys))
        inact.append(_parse_multiplicity(
            [x[1:-1] for x in elements if x.startswith('(') and x.endswith(')')],
            substance_keys
        ))

def as_RateExpr(self, unique_keys=None, constants=None, units=None, backend=math):
        from .rates import Eyring, MassAction
        args = [self.kB_h_times_exp_dS_R(constants, units, backend),
                self.dH_over_R(constants, units)]
        return MassAction(Eyring(args, unique_keys))

nargs = None
    parameter_keys = Radiolytic.parameter_keys + TPoly.parameter_keys

    def g_value(self, variables, backend):
        return self.eval_poly(variables, backend)


class RTPolyRadiolytic(RTPoly, Radiolytic):
    nargs = None
    parameter_keys = Radiolytic.parameter_keys + RTPoly.parameter_keys

    def g_value(self, variables, backend):
        return self.eval_poly(variables, backend)


class RTPolyMassAction(RTPoly, MassAction):
    """ Arguments: temperature_offset, c0, c1, ... """
    parameter_keys = RTPoly.parameter_keys
    nargs = None

    def rate_coeff(self, variables, backend=math):
        return self.eval_poly(variables, backend)


class _Log10XPolyMassAction(MassAction):
    skip_poly = 1  # kunit

    def rate_coeff(self, variables, backend=math):
        k_unit = self.arg(variables, 0)
        return 10**self.eval_poly(variables, backend)*k_unit

        return lambda *args, **kwargs: MassAction(Wrapper(*args, **kwargs))

def _r(r, p, substmap, parmap, *, unit_conc, unit_time,
       variables=None):
    pk, = r.param.unique_keys
    if isinstance(r.param, MassAction):
        ratcoeff = to_unitless(
            p[pk], unit_conc**(1-r.order())/unit_time
        )
        if not r.inact_reac:
            r_str = '{}, {}'.format(parmap[pk], r.string(substances=substmap, with_param=False,
                                                         Reaction_arrow='-->', Reaction_coeff_space=''))
        else:
            all_keys = r.keys()
            reac_stoichs = r.all_reac_stoich(all_keys)
            act_stoichs = r.active_reac_stoich(all_keys)
            rate = '*'.join([parmap[pk]] + [('%s^%d' % (substmap[k], v)) if v > 1 else substmap[k]
                                            for k, v in zip(all_keys, act_stoichs) if v > 0])
            r2 = Reaction(dict([(k, v) for k, v in zip(all_keys, reac_stoichs) if v]), r.prod)
            r_str = '{}, {}'.format(rate, r2.string(substances=substmap, with_param=False,
                                                    Reaction_arrow='\u21D2', Reaction_coeff_space=''))
    elif isinstance(r.param, RadiolyticBase):