function m = b2add(m, inFrom, stateFrom, 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(inFrom, 1, m.nU,  'kronecker:modelerror', 'Input Reactant #%d is out of range',    inFrom);  
    assertBetween(stateFrom, 1, m.nX,  'kronecker:modelerror', 'Reactant #%d is out of range',    stateFrom);  
    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 cooefecent is %d but must be +- 1,2', stoich);   
% Old way
%     [m.spinds.i_b2P m.spinds.j_b2P m.spinds.x_b2P m.spinds.n_b2P] = sparse_vec(m.spinds.i_b2P, m.spinds.j_b2P, m.spinds.x_b2P, m.spinds.n_b2P, (param-1)*m.nX + to, (inFrom-1)*m.nX+stateFrom, stoich);   
%     [m.spinds.i_b2 m.spinds.j_b2 m.spinds.x_b2 m.spinds.n_b2] = sparse_vec(m.spinds.i_b2, m.spinds.j_b2, m.spinds.x_b2, m.spinds.n_b2, to, (inFrom-1)*m.nX+stateFrom, stoich * m.p(param));


% Official way
    [m.spinds.i_b2P m.spinds.j_b2P m.spinds.x_b2P m.spinds.n_b2P] = sparse_vec(m.spinds.i_b2P, m.spinds.j_b2P, m.spinds.x_b2P, m.spinds.n_b2P, (param-1)*m.nX + to, (stateFrom-1)*m.nU+inFrom, stoich);   
    [m.spinds.i_b2 m.spinds.j_b2 m.spinds.x_b2 m.spinds.n_b2] = sparse_vec(m.spinds.i_b2, m.spinds.j_b2, m.spinds.x_b2, m.spinds.n_b2, to, (stateFrom-1)*m.nU+inFrom, stoich * m.p(param));

%     m.b2P((param-1)*m.nX + to, (inFrom-1)*m.nX+stateFrom) = ...
%         m.b2P((param-1)*m.nX + to, (inFrom-1)*m.nX+stateFrom) + stoich;        
%     m.b2(to, (inFrom-1)*m.nX+stateFrom) = ...
%         m.b2(to, (inFrom-1)*m.nX+stateFrom) + m.p(param) * stoich;
end