diff --git a/nextflow/scripts/discrete_symmetric.lphy b/nextflow/scripts/discrete_symmetric.lphy
new file mode 100644
index 0000000..5607b75
--- /dev/null
+++ b/nextflow/scripts/discrete_symmetric.lphy
@@ -0,0 +1,42 @@
+data {
+  // Specify options for reading the date information from file
+  options = {ageDirection="forward", ageRegex=".*_.*_(\d*\.\d+|\d+\.\d*)$"};
+  D = readNexus(file="data/simulated.nex", options=options);
+
+  // Retrieve number of sites in alignment and taxa names
+  L = D.nchar();
+  taxa = D.taxa();
+
+  // Extract Trait information and count # of traits
+  D_trait = extractTrait(taxa=taxa, sep="_", i=0);
+  K = D_trait.canonicalStateCount();
+
+  // Specify the data type for traits and number of pairwise combos
+  dim = K*(K-1)/2;
+  dataType = D_trait.dataType();
+}
+model {
+  // Model evolutionary history
+  Π ~ Dirichlet(conc=[2.0, 2.0, 2.0, 2.0]);      // Nucleotide Frequency prior
+  κ ~ LogNormal(meanlog=1.0, sdlog=1.25);        // Transition/transversion ratio prior
+  γ ~ LogNormal(meanlog=0.0, sdlog=2.0);         // Shape parameter for discretized gamma
+  r ~ DiscretizeGamma(shape=γ, ncat=4, 
+        replicates=L);                           // Site rates
+  Θ ~ LogNormal(meanlog=0.0, sdlog=1.0);         // Effective population size prior
+  ψ ~ Coalescent(taxa=taxa, theta=Θ);            // Coalescent time scaled phylogenetic tree prior
+  D ~ PhyloCTMC(Q=hky(kappa=κ, freq=Π), mu=0.004, 
+        siteRates=r, tree=ψ);                    // MCMC for evolution of alignment
+
+  // Model the geographic history
+  // modeled via P(t) = e^(Dt)
+  π_trait ~ Dirichlet(conc=rep(element=3.0, times=K));   // Base frequencies, assuming symmetry
+  I ~ Bernoulli(minSuccesses=dim-2, p=0.5, 
+        replicates=dim);                                 // Determines which rates are 0. 
+                                                         // This, plus select, implements BSSVS.
+  R_trait ~ Dirichlet(conc=rep(element=1.0, times=dim)); // Off-diagonal entries of the rate matrix.
+  Q_trait = generalTimeReversible(rates=select(x=R_trait, 
+              indicator=I), freq=π_trait);               // intantaneous rate matrix
+  μ_trait ~ LogNormal(meanlog=0, sdlog=1.25);            // migration events per unit time
+  D_trait ~ PhyloCTMC(L=1, Q=Q_trait, dataType=dataType, 
+              mu=μ_trait, tree=ψ);                       // MCMC for discrete locations
+}
\ No newline at end of file