function m = b1add(m, inFrom, 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(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);   
    
    [m.spinds.i_b1P m.spinds.j_b1P m.spinds.x_b1P m.spinds.n_b1P] = sparse_vec(m.spinds.i_b1P, m.spinds.j_b1P, m.spinds.x_b1P, m.spinds.n_b1P, (param-1)*m.nX + to, inFrom, stoich);
    
    [m.spinds.i_b1 m.spinds.j_b1 m.spinds.x_b1 m.spinds.n_b1] = sparse_vec(m.spinds.i_b1, m.spinds.j_b1, m.spinds.x_b1, m.spinds.n_b1, to, inFrom, stoich * m.p(param));

%     m.b1P((param-1)*m.nX + to,  inFrom) =  ...
%         m.b1P((param-1)*m.nX + to,  inFrom) + stoich;
% 
%     m.b1(to, inFrom) = m.b1(to, inFrom) + m.p(param)*stoich;
end