function m = a1add(m, from, to, param, stoich)
% 
% This work is licensed under the Creative Commons Attribution-Noncommercial-No Derivative 
% Works 3.0 United States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-nd/3.0/us/ or send a letter to Creative Commons,
% 171 Second Street, Suite 300, San Francisco, California, 94105, USA.
    assert(nargout == 1, 'kronecker:modelerror', 'function must be called with 1 output argument');
    assertBetween(from, 1, m.nX,  'kronecker:modelerror', 'Reactant #%d is out of range',    from);  
    assertBetween(to,   1, m.nX,  'kronecker:modelerror', 'Reactant #%d is out of range',      to);
    assertBetween(param, 1, m.nP, 'kronecker:modelerror', 'Parameter #%d is out of range',  param);
    assert(any(stoich == [-2, -1, 1, 2]), 'kronecker:modelerror', ...
        'Stoichiometric coefficient is %d but must be +- 1,2', stoich);

    [m.spinds.i_a1P m.spinds.j_a1P m.spinds.x_a1P m.spinds.n_a1P] = ...
        sparse_vec(m.spinds.i_a1P, m.spinds.j_a1P, m.spinds.x_a1P, m.spinds.n_a1P, (param-1)*m.nX + to, from, stoich);
    
    [m.spinds.i_a1 m.spinds.j_a1 m.spinds.x_a1 m.spinds.n_a1] = ...
        sparse_vec(m.spinds.i_a1, m.spinds.j_a1, m.spinds.x_a1, m.spinds.n_a1, to, from, stoich * m.p(param));
    
%     m.a1P((param-1)*m.nX + to,  from) =  m.a1P((param-1)*m.nX + to,  from) + stoich;
%     m.a1(to, from)                    =  m.a1(to, from) + stoich * m.p(param);
end