Files upload

Author: Huy Vo

DataSynthesizer/data_gen.m

@@ -0,0 +1,31 @@
+function datcell = data_gen( stoich_mat, propensity, x0, t_array, n_sample, unobservable)
+% datcell = data_gen( stoich_mat, propensity, x0, t_array, n_sample,
+% unobservable) synthesizes smFISH data at specific measurement time points
+% (for numerical testing purpose).
+datcell = cell( length(t_array ), 1 );
+n_species = length(x0);
+
+for it = 1:length( t_array )
+
+  t_ssa = t_array(it);
+
+
+  ssa_samples = zeros(n_sample, n_species);
+  for i = 1:n_sample
+    ssa_samples(i, :) = ssa(t_ssa, x0, stoich_mat, propensity);
+  end
+
+  ssa_samples = ssa_samples(:, setdiff((1:n_species), unobservable));
+
+  [observed_states, ~, i_observed] = unique(ssa_samples, 'rows');
+  n_observed = size(observed_states, 1);
+
+  [wbin, jbin] = hist(i_observed, unique(i_observed));
+
+  observed_freq = zeros(n_observed, 1);
+  observed_freq(jbin) = wbin;
+
+  datcell{it} = [observed_states observed_freq];
+end
+
+end

DataSynthesizer/ssa.m

@@ -0,0 +1,36 @@
+function end_state = ssa(t_f, x0, stoich_mat, propensity)
+% SSA stochastic simulation algorithm.
+%
+% t_f(in): final time.
+%
+% x0 (in): starting intial vector.
+%
+% stoich_mat (in): stoichiometry matrix of size (n.o.reactions) x
+% (n.o.species).
+%
+% propensity (in): function handle, propensity(x) = [a_1(x) a_2(x) .. a_M(x)]
+% where a_k(x) is the k-th propensity function evaluated at x.
+%
+% Reference:
+% =========
+% ?Gillespie, D. T. (1977). Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys. Chem., 81(25), 2340?2361. 
+%
+  x = x0;
+  t = 0;
+  while (t < t_f)
+    a = propensity(x);
+    a0 = sum(a);
+
+    u1 = rand();
+    u2 = rand();
+
+    % Update time
+    tau = min(log(1/u1)/a0, t_f - t);
+    t = t + tau;
+
+    % Update state
+    r = find(cumsum(a) > u2*a0, 1, 'first');
+    x = x + stoich_mat(r, :);
+  end
+  end_state = x;
+end

DataSynthesizer/test.m

@@ -0,0 +1,36 @@
+% Test the data generation on a toggle swtich model
+%% Define the model
+nmax = [50 50];
+
+x0 = [31 5];
+Pfull0 = zeros(prod(nmax+1),1); i0 = sub2ind(nmax+1, x0(1)+1, x0(2)+1); Pfull0(i0) = 1;
+
+S = [1 0;1 0;-1 0; 0 1; 0 1; 0 -1]'; % stoichiometry matrix
+
+ayx = 6.1e-3; axy = 2.6e-3; nyx = 2.1; nxy = 3; % we assume these parameters are known
+
+ind_prop = @(x,y) [ones(length(x),1), 1./(1+ayx*(y.^nyx)), x, ones(length(y),1), 1./(1+axy*(x.^nxy)), y]; % the parameter-independent part of the propensities
+
+params0 = [5e-5, 1e-2, 5e-4, 2e-3, 1e-2, 3e-4];
+
+t_array = linspace( 1, 10000, 4 );
+
+prop = @(x) params0.*ind_prop(x(1), x(2));
+datcell = data_gen( S', prop, x0, t_array, 1000, [] );
+%% Solve the FSP to compare with the generated data
+Acell = fsp_get_matrices( S, ind_prop, nmax );
+Pcell = fsp_solve( params0, t_array, Acell, Pfull0 );
+
+%% plot and compare results
+
+figure();
+for it = 1:length(t_array)
+    subplot( length(t_array), 2, 2*(it-1)+1);
+    contourf( reshape(Pcell{it}, nmax(1)+1, nmax(2)+1 ) );
+    subplot( length(t_array), 2, 2*(it-1)+2);
+    xdat = datcell{it}(:, 1:end-1);
+    wdat = datcell{it}(:,end);
+    contourf( full( sparse( xdat(:,1)+1, xdat(:,2)+1, wdat ) ) );
+    xlim ([0 nmax(2)+1] );
+    ylim ( [0 nmax(1)+1] );
+end