function [rxn_data, num_new, rxgp, cmpd_name, cmpd_nums, depend, min_depend, num_cmpds] = prepare_transcv2_rxn(rxn, rxgp, cmpd_name, cmpd_nums)

% need to pass in the number of copies of this gene.  Then do a loop for
% each copy with name mangling (rxn.gene = geneGSHA1, geneGSHA2 ... geneGSHA24)

% transcription reactions don't need to update probabilities
% unless the gene precursor has been made.  XTP change isn't enough, and would waste much cpu time.

% a non-zero 'min_depend' value for reaction A means that no other reactions will point into
% reaction A to cause an update of the probability of reaction A.
% However reactions created by this file can point into reaction A (and
% therefore cause an update of reaction A probability) via the 'depend' cell
% array.

% a zero value for min_depend for say reaction B indicates that the dependency matrix operates
% normally on reaction B -- create_dependency_matrix.m can find all other
% reactions in the reaction group that trigger an update of reaction B
% probability.


num_new = 0;
num_cmpds = length(cmpd_name);
num_bases = length(rxn.seq);
if rxn.polysize == 0
    rxn.polysize = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial binding and unbinding of polymerase to gene
% binding

num_new = num_new + 1;
rxn_data{num_new}.type = 'binding';
depend{num_new} = [3];
min_depend(num_new) = 0;

rxgp = build_comps(rxn.gene, rxgp, cmpd_name);
comp_id = find_comp_id(rxn.gene, rxgp.comps);
rxn_data{num_new}.substrates(1,1) = comp_id;

rxgp = build_comps(rxn.poly, rxgp, cmpd_name);
comp_id = find_comp_id(rxn.poly, rxgp.comps);
rxn_data{num_new}.substrates(2,1) = comp_id;

num_cmpds = num_cmpds + 1;
cmpd_name{num_cmpds} = [rxn.gene '_' rxn.poly '_closed'];
cmpd_nums(num_cmpds,:) = [0 0];
rxgp = build_comps(cmpd_name{num_cmpds}, rxgp, cmpd_name);
comp_id = find_comp_id(cmpd_name{num_cmpds},rxgp.comps);
rxn_data{num_new}.products(1,1) = comp_id;
rxn_data{num_new}.StoP.rate = rxn.rate(1);

% unbinding
if rxn.rate(2) ~= 0
    num_new = num_new + 1;
    rxn_data{num_new}.type = 'unbinding';
    depend{num_new} = [3];
    min_depend(num_new) = 0;

    rxn_data{num_new}.substrates(1,1) = comp_id;

    comp_id = find_comp_id(rxn.gene, rxgp.comps);
    rxn_data{num_new}.products(1,1) = comp_id;

    comp_id = find_comp_id(rxn.poly, rxgp.comps);
    rxn_data{num_new}.products(2,1) = comp_id;

    rxn_data{num_new}.StoP.rate = rxn.rate(2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% change closed to open form of gene_polymerase
sub1 = [rxn.gene '_o2'];
num_cmpds = num_cmpds + 1;
cmpd_name{num_cmpds} = sub1;
cmpd_nums(num_cmpds,:) = [1 0];

sub2 = [rxn.gene '_' rxn.poly '_closed'];

prod1 = [rxn.gene '_' rxn.poly '_open'];
num_cmpds = num_cmpds + 1;
cmpd_name{num_cmpds} = prod1;
cmpd_nums(num_cmpds,:) = [0 0];

num_new = num_new + 1;
if rxn.polysize == 1
    rxn_data{num_new}.type = 'unique-bi-bi';
    prod2 = rxn.gene;
    depend{num_new} = [num_new num_new+1];
else
    rxn_data{num_new}.type = 'binding';
    depend{num_new} = [num_new num_new+1];
end
min_depend(num_new) = 1;

rxgp = build_comps(sub1, rxgp, cmpd_name);
comp_id = find_comp_id(sub1, rxgp.comps);
rxn_data{num_new}.substrates(1,1) = comp_id;

%rxgp = build_comps(sub2, rxgp, cmpd_name);
comp_id = find_comp_id(sub2, rxgp.comps);
rxn_data{num_new}.substrates(2,1) = comp_id;

rxgp = build_comps(prod1, rxgp, cmpd_name);
comp_id = find_comp_id(prod1,rxgp.comps);
rxn_data{num_new}.products(1,1) = comp_id;

if rxn.polysize == 1
    rxgp = build_comps(prod2, rxgp, cmpd_name);
    comp_id = find_comp_id(prod2,rxgp.comps);
    rxn_data{num_new}.products(2,1) = comp_id;
end

rxn_data{num_new}.StoP.rate = rxn.rate(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bind 1st base
sub1 = [rxn.gene '_' rxn.poly '_open'];
if rxn.seq(1) == 'T'
    xtp = 'UTP';
else
    xtp = [rxn.seq(1) 'TP'];
end

prod1 = [rxn.gene '_1c'];
num_cmpds = num_cmpds + 1;
cmpd_name{num_cmpds} = prod1;
cmpd_nums(num_cmpds,:) = [0 0];

num_new = num_new + 1;
rxn_data{num_new}.type = 'binding';
depend{num_new} = [num_new num_new+1];
min_depend(num_new) = 1;

%rxgp = build_comps(sub1, rxgp, cmpd_name);
comp_id = find_comp_id(sub1, rxgp.comps);
rxn_data{num_new}.substrates(1,1) = comp_id;

rxgp = build_comps(xtp, rxgp, cmpd_name);
comp_id = find_comp_id(xtp, rxgp.comps);
rxn_data{num_new}.substrates(2,1) = comp_id;

rxgp = build_comps(prod1, rxgp, cmpd_name);
comp_id = find_comp_id(prod1,rxgp.comps);
rxn_data{num_new}.products(1,1) = comp_id;
rxn_data{num_new}.StoP.rate = rxn.rate(4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now the big loop for adding each base in 3 steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ig = 2:num_bases
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %bind second site
    sub1 = [rxn.gene '_' int2str(ig-1) 'c'];
    if rxn.seq(ig) == 'T'
        xtp = 'UTP';
    else
        xtp = [rxn.seq(ig) 'TP'];
    end
    prod1 = [rxn.gene '_' int2str(ig) 'a'];
    num_cmpds = num_cmpds + 1;
    cmpd_name{num_cmpds} = prod1;
    cmpd_nums(num_cmpds,:) = [0 0];
    
    num_new = num_new + 1;
    rxn_data{num_new}.type = 'binding';
    depend{num_new} = [num_new num_new+1];
    min_depend(num_new) = 1;

    rxgp = build_comps(sub1, rxgp, cmpd_name);
    comp_id = find_comp_id(sub1, rxgp.comps);
    rxn_data{num_new}.substrates(1,1) = comp_id;

    rxgp = build_comps(xtp, rxgp, cmpd_name);
    comp_id = find_comp_id(xtp, rxgp.comps);
    rxn_data{num_new}.substrates(2,1) = comp_id;

    rxgp = build_comps(prod1, rxgp, cmpd_name);
    comp_id = find_comp_id(prod1,rxgp.comps);
    rxn_data{num_new}.products(1,1) = comp_id;
    rxn_data{num_new}.StoP.rate = rxn.rate(5);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % link site 1 and site 2
    sub1 = [rxn.gene '_' int2str(ig) 'a'];
    prod1 = [rxn.gene '_' int2str(ig) 'b'];
    num_cmpds = num_cmpds + 1;
    cmpd_name{num_cmpds} = prod1;
    cmpd_nums(num_cmpds,:) = [0 0];
    prod2 = 'PPi';

    num_new = num_new + 1;
    rxn_data{num_new}.type = 'unbinding';
    depend{num_new} = [num_new num_new+1];
    min_depend(num_new) = 1;
    
    rxgp = build_comps(sub1, rxgp, cmpd_name);
    comp_id = find_comp_id(sub1, rxgp.comps);
    rxn_data{num_new}.substrates(1,1) = comp_id;

    rxgp = build_comps(prod1, rxgp, cmpd_name);
    comp_id = find_comp_id(prod1,rxgp.comps);
    rxn_data{num_new}.products(1,1) = comp_id;

    rxgp = build_comps(prod2, rxgp, cmpd_name);
    comp_id = find_comp_id(prod2,rxgp.comps);
    rxn_data{num_new}.products(2,1) = comp_id;
    rxn_data{num_new}.StoP.rate = rxn.rate(6);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % slide base at site 2 to site 1
    sub1 = [rxn.gene '_o' int2str(ig+1)];
    num_cmpds = num_cmpds + 1;
    cmpd_name{num_cmpds} = sub1;
    cmpd_nums(num_cmpds,:) = [1 0];
    sub2 = [rxn.gene '_' int2str(ig) 'b'];
    prod1 = [rxn.gene '_' int2str(ig) 'c'];
    num_cmpds = num_cmpds + 1;
    cmpd_name{num_cmpds} = prod1;
    cmpd_nums(num_cmpds,:) = [0 0];
    
    num_new = num_new + 1;
    if ig > rxn.polysize + 1
        rxn_data{num_new}.type = 'unique-bi-bi';
        prod2 = [rxn.gene '_o' int2str(ig - rxn.polysize + 1)];
        depend{num_new} = [num_new - (rxn.polysize.*3) num_new num_new+1];
    elseif ig == rxn.polysize + 1
        rxn_data{num_new}.type = 'unique-bi-bi';
        prod2 = [rxn.gene '_o' int2str(ig - rxn.polysize + 1)];
        depend{num_new} = [num_new - (rxn.polysize.*3) - 1 num_new num_new+1];
    elseif ig == rxn.polysize
        rxn_data{num_new}.type = 'unique-bi-bi';
        prod2 = [rxn.gene];
        depend{num_new} = [num_new num_new+1];
    else
        rxn_data{num_new}.type = 'binding';
        depend{num_new} = [num_new num_new+1];
    end
    min_depend(num_new) = 1;
    
    rxgp = build_comps(sub1, rxgp, cmpd_name);
    comp_id = find_comp_id(sub1, rxgp.comps);
    rxn_data{num_new}.substrates(1,1) = comp_id;

    rxgp = build_comps(sub2, rxgp, cmpd_name);
    comp_id = find_comp_id(sub2, rxgp.comps);
    rxn_data{num_new}.substrates(2,1) = comp_id;

    rxgp = build_comps(prod1, rxgp, cmpd_name);
    comp_id = find_comp_id(prod1,rxgp.comps);
    rxn_data{num_new}.products(1,1) = comp_id;
    
    if ig >= rxn.polysize
        rxgp = build_comps(prod2, rxgp, cmpd_name);
        comp_id = find_comp_id(prod2,rxgp.comps);
        rxn_data{num_new}.products(2,1) = comp_id;
    end
    
    rxn_data{num_new}.StoP.rate = rxn.rate(7);
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% release mRNA, release polymerase, release all the sites that were blocked
% by the polymerase


sub1 = [rxn.gene '_' int2str(num_bases) 'c'];
prod1 = rxn.mRNA;
prod2 = rxn.poly;

num_new = num_new + 1;
rxn_data{num_new}.type = 'uni-multi';
depend{num_new} = [num_new];
min_depend(num_new) = 1;

rxgp = build_comps(sub1, rxgp, cmpd_name);
comp_id = find_comp_id(sub1, rxgp.comps);
rxn_data{num_new}.substrates(1,1) = comp_id;

rxgp = build_comps(prod1, rxgp, cmpd_name);
comp_id = find_comp_id(prod1,rxgp.comps);
rxn_data{num_new}.products(1,1) = comp_id;
rxn_data{num_new}.products(1,2) = 1;
    
rxgp = build_comps(prod2, rxgp, cmpd_name);
comp_id = find_comp_id(prod2,rxgp.comps);
rxn_data{num_new}.products(2,1) = comp_id;
rxn_data{num_new}.products(2,2) = 1;

other_prod_num = 2;
for ig =  num_bases + 1 : -1 : (num_bases - rxn.polysize + 2)
    other_prod_num = other_prod_num + 1;
    if ig ~= 1
        prodn = [rxn.gene '_o' int2str(ig)];
    else
        prodn = rxn.gene;
    end
    rxgp = build_comps(prodn, rxgp, cmpd_name);
    comp_id = find_comp_id(prodn,rxgp.comps);
    rxn_data{num_new}.products(other_prod_num,1) = comp_id;
    rxn_data{num_new}.products(other_prod_num,2) = 1;
    
    if ig >= 3
        depend{num_new} = [num_new - (num_bases+1-ig)*3 - 1 depend{num_new}];
    elseif ig == 2
        depend{num_new} = [num_new - (num_bases+1-ig)*3 - 2 depend{num_new}];
    elseif ig == 1
        break;
    end
    
    
end
rxn_data{num_new}.StoP.rate = rxn.rate(8);


