function m = symbolic2MassAction(SymModel, uNames, yNames, yMembers, opts)

verbose = logical(opts.Verbose);
opts.Verbose = max(opts.Verbose-1, 0);

name        = SymModel.Name;
xNames      = SymModel.Xnames;
kNames      = SymModel.Knames;
x0          = SymModel.x0;
k           = SymModel.k;
xSyms       = SymModel.Xsyms;
kSyms       = SymModel.Ksyms;
rSyms       = SymModel.Rsyms;
S           = SymModel.S;

%% Parse the input
% Standardize uNames as a cell array of strings
if isa(uNames, 'char')
    uNames = {uNames};
    Uind = find(strcmp(xNames, uNames{1}));
elseif isa(uNames, 'sym')
    temp = uNames;
    uNames = cell(length(uNames),1);
    for i = 1:length(uNames)
        uNames{i} = char(temp(i));
        Uind = find(strcmp(xNames, uNames{i}));
    end
elseif isa(uNames, 'numeric')
    Uind = uNames;
    uNames = xNames(Uind);
elseif isa(uNames, 'cell')
    Uind = zeros(length(uNames),1);
    for i = 1:length(uNames)
        Uind(i) = find(strcmp(xNames, uNames{i}));
    end
end

% Remove inputs from x
uSyms           = xSyms(Uind);
xNames(Uind)    = [];
xSyms(Uind)     = [];
x0(Uind)        = [];
S(Uind,:)       = [];

%% Basic constants
% Basic constants
nu = numel(uNames);
nx = numel(xNames);
nk = numel(kNames);

%% Parse the output
% Standardize  yNames as a cell array of strings
% and yMembers as a cell array of cell arrays of strings
if isempty(yMembers)
    % yNames are the yMembers
    if isempty(yNames)
        % All species are outputs
        yNames = 1:nx;
    end
    
    if isa(yNames, 'char')
        ny = 1;
        yNames = {yNames};
        C = sparse(ny, nx);
        C(strcmp(xNames, yNames{1})) = 1;
    elseif isa(yNames, 'sym')
        ny = length(yNames);
        temp = yNames;
        yNames = cell(ny,1);
        C = sparse(ny, nx);
        for i = 1:ny
            yNames{i} = char(temp(i));
            C(i, strcmp(xNames, yNames{i})) = C(i, strcmp(xNames, yNames{i})) + 1;
        end
    elseif isa(yNames, 'numeric')
        ny = length(yNames);
        yNames = xNames(yNames);
        C = sparse(ny, nx);
        for i = 1:ny
            C(i, strcmp(xNames, yNames{i})) = C(i, strcmp(xNames, yNames{i})) + 1;
        end
    else%if isa(yNames, 'cell')
        ny = length(yNames);
        yNames = vec(yNames);
        C = sparse(ny, nx);
        for i = 1:ny
            C(i, strcmp(xNames, yNames{i})) = C(i, strcmp(xNames, yNames{i})) + 1;
        end
    end
else
    % yMembers is provided
    if isa(yNames, 'char')
        ny = 1;
        yNames = {yNames};
    elseif isa(yNames, 'sym')
        ny = length(yNames);
        temp = yNames;
        yNames = cell(ny,1);
        for i = 1:ny
            yNames{i} = char(temp(i));
        end
    elseif isa(yNames, 'numeric')
        ny = length(yNames);
        yNames = xNames(yNames);
    else%if isa(yNames, 'cell')
        ny = length(yNames);
        yNames = vec(yNames);
    end
    
    % yMembers must be a cell array of cells
    C = sparse(ny, nx);
    for i = 1:ny
        nMem = length(yMembers{i});
        for j = 1:nMem
            C(i, strcmp(xNames, yMembers{i}{j})) = C(i, strcmp(xNames, yMembers{i}{j})) + 1;
        end
    end
end

%% Sanity check
assert(all(vec(fix(S) == S)) && all(vec(abs(S) <= 2)), 'KroneckerBio:symbolic2MassAction:InvalidStoichiometry', 'Invalid stoichiometry matrix does not consist entirely of -2, -1, 0, 1, and 2')

%% Initialize model
% Basic structure
m = initModel(nx, nu, ny, nk);

m.name = name;
m.xNames = xNames;
m.uNames = uNames;
m.yNames = yNames;
m.pNames = kNames;
m.p = k;
m.ic = x0;
m.c = C;

%% Loop over reaction rate equations
if verbose; fprintf('Parsing reaction symbolics to kronecker reactions...\n'); end
for iEqu = 1:numel(rSyms)
    % Extract all variables used in this reaction rate
    allvar = symvar(rSyms(iEqu));
    
    % Variables that are rate constants; index in kSyms
    kPosEqu = lookup(allvar, kSyms);
    kPosEqu = kPosEqu(kPosEqu ~= 0);
    
    % Number of microscopic reactions in this rate
    nMic = numel(kPosEqu);
    
    % Count of products and reactants in stoichiometry matrix
    % Double up on species that get twice as much
    reac = [find(S(:,iEqu) < 0); find(S(:,iEqu) == -2)];
    prod = [find(S(:,iEqu) > 0); find(S(:,iEqu) == 2)];
    nreac = numel(reac);
    nprod = numel(prod);
    assert(nreac <= 2 && nprod <= 2, 'KroneckerBio:symbolic2MassAction:InvalidStoichiometry', 'Reaction rate %i has an unacceptable stoichiometry vector', iEqu)
    
    % Loop over each microscopic reaction and add it
    for iMic = 1:nMic
        % Differentiate wrt the rate constant to return only the species
        derMic = diff(rSyms(iEqu), kSyms(kPosEqu(iMic)));
        
        % Extract those species involved
        allx = symvar(derMic);
        
        % Ensure second order or less with inputs and species
        xuPosMic = lookup(allx, [xSyms; uSyms]);
        xPosMic = xuPosMic(xuPosMic <= nx);
        uPosMic = xuPosMic(xuPosMic > nx) - nx;
        nxuMic = numel(xuPosMic);
        nxMic = numel(xPosMic);
        nuMic = numel(uPosMic);
        assert(all(xuPosMic ~= 0) && nxuMic <= 2, 'KroneckerBio:symbolic2MassAction:InvalidReaction', 'Reaction rate %i does not represent a quadratic ODE', iEqu)
        
        % Switch on input
        if nuMic
            % Reaction has an input
            assert(nuMic < 2, 'KroneckerBio:symbolic2MassAction:SecondOrderInput', 'Kronecker Bio does not currently support second-order reactions between inputs (see Reaction %i); create a new input equal to the product of those two inputs', iEqu)
            
            % Switch on the number of states involved
            switch(nxMic)
                case(0)
                    % Derivative check for first-order reaction
                    derMic = diff(derMic, uSyms(uPosMic));
                    assert(derMic == -1 || derMic == 1, 'KroneckerBio:symbolic2MassAction:ZerothOrderCheck', 'Reaction rate %i has a zeroth-order term that failed a sanity check', iEqu)
                    
                    if derMic == 1
                        % Add first-order input forward
                        assert(nreac == 0, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a first-order input term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of products
                        switch(nprod)
                            case(0)
                                % Pointless reaction
                                % m = addInputOne2Zero(m, kPosEqu(iMic));
                            case(1)
                                m = addInputOne2One(m, uPosMic, prod, kPosEqu(iMic));
                            case(2)
                                m = addInputOne2Two(m, uPosMic, min(prod), max(prod), kPosEqu(iMic));
                        end
                    elseif derMic == -1
                        % Add first-order input reverse
                        assert(nprod == 0, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a first-order input term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of "reactants"
                        switch(nreac)
                            case(0)
                                % Pointless reaction
                                %m = addInputOne2Zero(m, uPosMic, kPosEqu(iMic));
                            case(1)
                                m = addInputOne2One(m, uPosMic, reac, kPosEqu(iMic));
                            case(2)
                                m = addInputOne2Two(m, uPosMic, min(reac), max(reac), kPosEqu(iMic));
                        end
                    end
                case(1)
                    % Derivative check for first-order reaction
                    derMic = diff(derMic, uSyms(uPosMic));
                    derMic = diff(derMic, xSyms(xPosMic));
                    assert(derMic == -1 || derMic == 1, 'KroneckerBio:symbolic2MassAction:ZerothOrderCheck', 'Reaction rate %i has a second-order input term that failed a sanity check', iEqu)
                    
                    if derMic == 1
                        % Add second-order input forward
                        assert(nreac == 1, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a second-order input term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of products
                        switch(nprod)
                            case(0)
                                m = addInputTwo2Zero(m, uPosMic, xPosMic, kPosEqu(iMic));
                            case(1)
                                m = addInputTwo2One(m, uPosMic, xPosMic, prod, kPosEqu(iMic));
                            case(2)
                                m = addInputTwo2Two(m, uPosMic, xPosMic, min(prod), max(prod), kPosEqu(iMic));
                        end
                    elseif derMic == -1
                        % Add second-order input reverse
                        assert(nprod == 1, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a second-order input term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of "reactants"
                        switch(nreac)
                            case(0)
                                m = addInputTwo2Zero(m, uPosMic, xPosMic, kPosEqu(iMic));
                            case(1)
                                m = addInputTwo2One(m, uPosMic, xPosMic, reac, kPosEqu(iMic));
                            case(2)
                                m = addInputTwo2Two(m, uPosMic, xPosMic, min(reac), max(reac), kPosEqu(iMic));
                        end
                    end
            end
        else
            % Reaction has no input
            % Switch one zeroth, first, or second order
            switch(nxMic)
                case(0)
                    % Zeroth-order reaction
                    assert(derMic == -1 || derMic == 1, 'KroneckerBio:symbolic2MassAction:ZerothOrderCheck', 'Reaction rate %i has a zeroth-order term that failed a sanity check', iEqu)
                    
                    if derMic == 1
                        % Add zeroth order synthesis forward
                        assert(nreac == 0, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a zeroth-order term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of products
                        switch(nprod)
                            case(0)
                                % Pointless reaction
                                % m = addZero2Zero(m, kPosEqu(iMic));
                            case(1)
                                m = addZero2One(m, prod, kPosEqu(iMic));
                            case(2)
                                m = addZero2Two(m, min(prod), max(prod), kPosEqu(iMic));
                        end
                    elseif derMic == -1
                        % Add zeroth order synthesis reverse
                        assert(nprod == 0, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a zeroth-order term that does not match the stoichiometry', iEqu)

                        % Switch on number of "reactants"
                        switch(nreac)
                            case(0)
                                % Pointless reaction
                                % m = addZero2Zero(m, kPosEqu(iMic));
                            case(1)
                                m = addZero2One(m, reac, kPosEqu(iMic));
                            case(2)
                                m = addZero2Two(m, min(reac), max(reac), kPosEqu(iMic));
                        end
                    end
                case(1)
                    % Derivative check for first-order reaction
                    derMic = diff(derMic, xSyms(xPosMic));
                    
                    % Switch on actual first order or self second-order
                    if derMic == -1 || derMic == 1
                        % True first-order reaction
                        
                        if derMic == 1
                            % Add first order forward
                            assert(nreac == 1, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a single species term that does not match the stoichiometry', iEqu)
                            
                            % Switch on number of products
                            switch(nprod)
                                case(0)
                                    m = addOne2Zero(m, xPosMic, kPosEqu(iMic));
                                case(1)
                                    m = addOne2One(m, xPosMic, prod, kPosEqu(iMic));
                                case(2)
                                    m = addOne2Two(m, xPosMic, min(prod), max(prod), kPosEqu(iMic));
                            end
                        elseif derMic == -1
                            % Add first order reverse
                            assert(nprod == 1, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a single species term that does not match the stoichiometry', iEqu)

                            % Switch on number of "reactants"
                            switch(nreac)
                                case(0)
                                    m = addOne2Zero(m, xPosMic, kPosEqu(iMic));
                                case(1)
                                    m = addOne2One(m, xPosMic, reac, kPosEqu(iMic));
                                case(2)
                                    m = addOne2Two(m, xPosMic, min(reac), max(reac), kPosEqu(iMic));
                            end
                        end
                    else
                        % Species reacts with itself
                        derMic = diff(derMic, xSyms(xPosMic));
                        assert(derMic == 2 || derMic == -2, 'KroneckerBio:symbolic2MassAction:FirstOrderCheck', 'Reaction rate %i has a single species term that failed a sanity check', iEqu)

                        if derMic == 2
                            % Add second order forward
                            assert(nreac == 2, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a single species term that does not match the stoichiometry', iEqu)
                            
                            % Switch on number of products
                            switch(nprod)
                                case(0)
                                    m = addTwo2Zero(m, xPosMic, xPosMic, kPosEqu(iMic));
                                case(1)
                                    m = addTwo2One(m, xPosMic, xPosMic, prod, kPosEqu(iMic));
                                case(2)
                                    m = addTwo2Two(m, xPosMic, xPosMic, min(prod), max(prod), kPosEqu(iMic));
                            end
                        elseif derMic == -2
                            % Add second order reverse
                            assert(nprod == 2, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a single species term that does not match the stoichiometry', iEqu)
                            
                            % Switch on number of "reactants"
                            switch(nreac)
                                case(0)
                                    m = addTwo2Zero(m, xPosMic, xPosMic, kPosEqu(iMic));
                                case(1)
                                    m = addTwo2One(m, xPosMic, xPosMic, reac, kPosEqu(iMic));
                                case(2)
                                    m = addTwo2Two(m, xPosMic, xPosMic, min(reac), max(reac), kPosEqu(iMic));
                            end
                        end
                    end
                case(2)
                    % Derivative check for second-order reaction
                    derMic = diff(derMic, xSyms(xPosMic(1)));
                    derMic = diff(derMic, xSyms(xPosMic(2)));
                    assert(derMic == 1 || derMic == -1, 'KroneckerBio:symbolic2MassAction:SecondOrderCheck', 'Reaction rate %i has a double species term that failed a sanity check', iEqu)
                    
                    if derMic == 1
                        % Add second order forward
                        assert(nreac == 2, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a double species term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of products
                        switch(nprod)
                            case(0)
                                m = addTwo2Zero(m, xPosMic(1), xPosMic(2), kPosEqu(iMic));
                            case(1)
                                m = addTwo2One(m, xPosMic(1), xPosMic(2), prod, kPosEqu(iMic));
                            case(2)
                                m = addTwo2Two(m, xPosMic(1), xPosMic(2), min(prod), max(prod), kPosEqu(iMic));
                        end
                    elseif derMic == -1
                        % Add second order reverse
                        assert(nprod == 2, 'KroneckerBio:symbolic2MassAction:StoichiometryMismatch', 'Reaction rate %i has a double species term that does not match the stoichiometry', iEqu)
                        
                        % Switch on number of "reactants"
                        switch(nreac)
                            case(0)
                                m = addTwo2Zero(m, xPosMic(1), xPosMic(2), kPosEqu(iMic));
                            case(1)
                                m = addTwo2One(m, xPosMic(1), xPosMic(2), reac, kPosEqu(iMic));
                            case(2)
                                m = addTwo2Two(m, xPosMic(1), xPosMic(2), min(reac), max(reac), kPosEqu(iMic));
                        end
                    end
                    
            end
        end
    end
end

%% Finalize model
m = postProcessModel(m);
