function m = a2add(m,from1, from2, 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(from1, 1, m.nX,  'kronecker:modelerror', 'Reactant #%d is out of range',    from1);  
    assertBetween(from2, 1, m.nX,  'kronecker:modelerror', 'Reactant #%d is out of range',    from2);  
    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_a2P m.spinds.j_a2P m.spinds.x_a2P m.spinds.n_a2P] = ...
        sparse_vec(m.spinds.i_a2P, m.spinds.j_a2P, m.spinds.x_a2P, m.spinds.n_a2P, (param-1)*m.nX + to, (from1-1)*m.nX+from2, stoich/2);
    [m.spinds.i_a2P m.spinds.j_a2P m.spinds.x_a2P m.spinds.n_a2P] = ...
        sparse_vec(m.spinds.i_a2P, m.spinds.j_a2P, m.spinds.x_a2P, m.spinds.n_a2P, (param-1)*m.nX + to, (from2-1)*m.nX+from1, stoich/2);
    
    [m.spinds.i_a2 m.spinds.j_a2 m.spinds.x_a2 m.spinds.n_a2] = ...
        sparse_vec(m.spinds.i_a2, m.spinds.j_a2, m.spinds.x_a2, m.spinds.n_a2, to, (from1-1)*m.nX+from2, stoich * m.p(param)/2);
    [m.spinds.i_a2 m.spinds.j_a2 m.spinds.x_a2 m.spinds.n_a2] = ...
        sparse_vec(m.spinds.i_a2, m.spinds.j_a2, m.spinds.x_a2, m.spinds.n_a2, to, (from2-1)*m.nX+from1, stoich * m.p(param)/2);
    
%     % Enter the stoich symmetrically into the matrix
%     m.a2P((param-1)*m.nX + to, (from1-1)*m.nX+from2) = ...
%     m.a2P((param-1)*m.nX + to, (from1-1)*m.nX+from2) + stoich/2;
% 
%     m.a2P((param-1)*m.nX + to, (from2-1)*m.nX+from1) = ...
%     m.a2P((param-1)*m.nX + to, (from2-1)*m.nX+from1) + stoich/2;
% 
%     % Enter the rate symmetrically into the matrix
%     m.a2(to, (from1-1)*m.nX+from2) = ...
%     m.a2(to, (from1-1)*m.nX+from2) + stoich * m.p(param)/2;
% 
%     m.a2(to, (from2-1)*m.nX+from1) = ...
%     m.a2(to, (from2-1)*m.nX+from1) + stoich * m.p(param)/2;
end