diff --git a/docs/build/index.html b/docs/build/index.html
index 991dc9f..b340372 100644
--- a/docs/build/index.html
+++ b/docs/build/index.html
@@ -8,7 +8,7 @@
 @push acs 1.0 X ⟶ Y 
 @prob_init acs X=1 Y=2 XY=α
 @prob_params acs γ=1 α=4
-@solve_and_plot acs</code></pre></div></section></article><h2 id="Modify-a-model"><a class="docs-heading-anchor" href="#Modify-a-model">Modify a model</a><a id="Modify-a-model-1"></a><a class="docs-heading-anchor-permalink" href="#Modify-a-model" title="Permalink"></a></h2><p>We list common transition attributes:</p><table><tr><th style="text-align: left">attribute</th><th style="text-align: left">interpretation</th></tr><tr><td style="text-align: left"><code>transPriority</code></td><td style="text-align: left">priority of a transition (influences resource allocation)</td></tr><tr><td style="text-align: left"><code>transProbOfSuccess</code></td><td style="text-align: left">probability that a transition terminates successfully</td></tr><tr><td style="text-align: left"><code>transCapacity</code></td><td style="text-align: left">maximum number of concurrent instances of the transition</td></tr><tr><td style="text-align: left"><code>transCycleTime</code></td><td style="text-align: left">duration of a transition&#39;s instance (adjusted by resource allocation)</td></tr><tr><td style="text-align: left"><code>transMaxLifeTime</code></td><td style="text-align: left">maximal duration of a transition&#39;s instance</td></tr><tr><td style="text-align: left"><code>transPostAction</code></td><td style="text-align: left">action to be executed once a transition&#39;s instance terminates</td></tr><tr><td style="text-align: left"><code>transName</code></td><td style="text-align: left">name of a transition</td></tr></table><p>We list common species attributes:</p><table><tr><th style="text-align: left">attribute</th><th style="text-align: left">interpretation</th></tr><tr><td style="text-align: left"><code>specInitUncertainty</code></td><td style="text-align: left">uncertainty about variable&#39;s initial state (modelled as Gaussian standard deviation)</td></tr><tr><td style="text-align: left"><code>specInitVal</code></td><td style="text-align: left">initial value of a variable</td></tr></table><p>Moreover, it is possible to specify the semantics of the &quot;rate&quot; term. By default, at each time step <code>n ~ Poisson(rate * dt)</code> instances of a given transition will be spawned. If you want to specify the rate in terms of a cycle time, you may want to use <code>@ct(cycle_time)</code>, e.g., <code>@ct(ex), A --&gt; B, ...</code>. This is a shorthand for <code>1/ex, A --&gt; B, ...</code>.</p><p>For deterministic &quot;rates&quot;, use <code>@per_step(ex)</code>. Here, <code>ex</code> evaluates to a deterministic number (ceiled to the nearest integer) of a transition&#39;s instances to spawn per a single integrator&#39;s step. However, note that in this case, the number doesn&#39;t scale with the step length! Moreover</p><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@add_species" href="#ReactiveDynamics.@add_species"><code>ReactiveDynamics.@add_species</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Add new species to a model.</p><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@add_species acs S I R</code></pre></div></section></article><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@aka" href="#ReactiveDynamics.@aka"><code>ReactiveDynamics.@aka</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Alias object name in an acs.</p><p><strong>Default names</strong></p><table><tr><th style="text-align: left">name</th><th style="text-align: left">short name</th></tr><tr><td style="text-align: left">species</td><td style="text-align: left">S</td></tr><tr><td style="text-align: left">transition</td><td style="text-align: left">T</td></tr><tr><td style="text-align: left">action</td><td style="text-align: left">A</td></tr><tr><td style="text-align: left">event</td><td style="text-align: left">E</td></tr><tr><td style="text-align: left">param</td><td style="text-align: left">P</td></tr><tr><td style="text-align: left">meta</td><td style="text-align: left">M</td></tr></table><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@aka acs species=resource transition=reaction</code></pre></div></section></article><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@mode" href="#ReactiveDynamics.@mode"><code>ReactiveDynamics.@mode</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Set species modality.</p><p><strong>Supported modalities</strong></p><ul><li>nonblock</li><li>conserved</li><li>rate</li></ul><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@mode acs (r&quot;proj\w+&quot;, r&quot;experimental\w+&quot;) conserved
+@solve_and_plot acs</code></pre></div></section></article><h2 id="Modify-a-model"><a class="docs-heading-anchor" href="#Modify-a-model">Modify a model</a><a id="Modify-a-model-1"></a><a class="docs-heading-anchor-permalink" href="#Modify-a-model" title="Permalink"></a></h2><p>We list common transition attributes:</p><table><tr><th style="text-align: left">attribute</th><th style="text-align: left">interpretation</th></tr><tr><td style="text-align: left"><code>transPriority</code></td><td style="text-align: left">priority of a transition (influences resource allocation)</td></tr><tr><td style="text-align: left"><code>transProbOfSuccess</code></td><td style="text-align: left">probability that a transition terminates successfully</td></tr><tr><td style="text-align: left"><code>transCapacity</code></td><td style="text-align: left">maximum number of concurrent instances of the transition</td></tr><tr><td style="text-align: left"><code>transCycleTime</code></td><td style="text-align: left">duration of a transition&#39;s instance (adjusted by resource allocation)</td></tr><tr><td style="text-align: left"><code>transMaxLifeTime</code></td><td style="text-align: left">maximal duration of a transition&#39;s instance</td></tr><tr><td style="text-align: left"><code>transPostAction</code></td><td style="text-align: left">action to be executed once a transition&#39;s instance terminates</td></tr><tr><td style="text-align: left"><code>transName</code></td><td style="text-align: left">name of a transition</td></tr></table><p>We list common species attributes:</p><table><tr><th style="text-align: left">attribute</th><th style="text-align: left">interpretation</th></tr><tr><td style="text-align: left"><code>specInitUncertainty</code></td><td style="text-align: left">uncertainty about variable&#39;s initial state (modelled as Gaussian standard deviation)</td></tr><tr><td style="text-align: left"><code>specInitVal</code></td><td style="text-align: left">initial value of a variable</td></tr></table><p>Moreover, it is possible to specify the semantics of the &quot;rate&quot; term. By default, at each time step <code>n ~ Poisson(rate * dt)</code> instances of a given transition will be spawned. If you want to specify the rate in terms of a cycle time, you may want to use <code>@ct(cycle_time)</code>, e.g., <code>@ct(ex), A --&gt; B, ...</code>. This is a shorthand for <code>1/ex, A --&gt; B, ...</code>.</p><p>For deterministic &quot;rates&quot;, use <code>@deterministic(ex)</code>. Here, <code>ex</code> evaluates to a deterministic number (ceiled to the nearest integer) of a transition&#39;s instances to spawn per a single integrator&#39;s step. However, note that in this case, the number doesn&#39;t scale with the step length! Moreover</p><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@add_species" href="#ReactiveDynamics.@add_species"><code>ReactiveDynamics.@add_species</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Add new species to a model.</p><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@add_species acs S I R</code></pre></div></section></article><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@aka" href="#ReactiveDynamics.@aka"><code>ReactiveDynamics.@aka</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Alias object name in an acs.</p><p><strong>Default names</strong></p><table><tr><th style="text-align: left">name</th><th style="text-align: left">short name</th></tr><tr><td style="text-align: left">species</td><td style="text-align: left">S</td></tr><tr><td style="text-align: left">transition</td><td style="text-align: left">T</td></tr><tr><td style="text-align: left">action</td><td style="text-align: left">A</td></tr><tr><td style="text-align: left">event</td><td style="text-align: left">E</td></tr><tr><td style="text-align: left">param</td><td style="text-align: left">P</td></tr><tr><td style="text-align: left">meta</td><td style="text-align: left">M</td></tr></table><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@aka acs species=resource transition=reaction</code></pre></div></section></article><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@mode" href="#ReactiveDynamics.@mode"><code>ReactiveDynamics.@mode</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Set species modality.</p><p><strong>Supported modalities</strong></p><ul><li>nonblock</li><li>conserved</li><li>rate</li></ul><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@mode acs (r&quot;proj\w+&quot;, r&quot;experimental\w+&quot;) conserved
 @mode acs (S, I) conserved
 @mode acs S conserved</code></pre></div></section></article><article class="docstring"><header><a class="docstring-binding" id="ReactiveDynamics.@name_transition" href="#ReactiveDynamics.@name_transition"><code>ReactiveDynamics.@name_transition</code></a> — <span class="docstring-category">Macro</span></header><section><div><p>Set name of a transition in the model.</p><p><strong>Examples</strong></p><pre><code class="language-julia hljs">@name_transition acs 1=&quot;name&quot;
 @name_transition acs name=&quot;transition_name&quot;
diff --git a/docs/build/search_index.js b/docs/build/search_index.js
index 6b6953e..7c23342 100644
--- a/docs/build/search_index.js
+++ b/docs/build/search_index.js
@@ -1,3 +1,3 @@
 var documenterSearchIndex = {"docs":
-[{"location":"index.html#API-Documentation","page":"API Documentation","title":"API Documentation","text":"","category":"section"},{"location":"index.html#Create-a-model","page":"API Documentation","title":"Create a model","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@ReactionNetwork","category":"page"},{"location":"index.html#ReactiveDynamics.@ReactionNetwork","page":"API Documentation","title":"ReactiveDynamics.@ReactionNetwork","text":"Macro that takes an expression corresponding to a reaction network and outputs an instance of TheoryReactionNetwork that can be converted to a DiscreteProblem or solved directly.\n\nMost arrows accepted (both right, left, and bi-drectional arrows). Use 0 or ∅ for annihilation/creation to/from nothing.\n\nCustom functions and sampleable objects can be used as numeric parameters. Note that these have to be accessible from ReactiveDynamics's source code.\n\nExamples\n\nacs = @ReactionNetworkSchema begin\n    1.0, X ⟶ Y\n    1.0, X ⟶ Y, priority=>6., prob=>.7, capacity=>3.\n    1.0, ∅ --> (Poisson(.3γ)X, Poisson(.5)Y)\n    (XY > 100) && (XY -= 1)\nend\n@push acs 1.0 X ⟶ Y \n@prob_init acs X=1 Y=2 XY=α\n@prob_params acs γ=1 α=4\n@solve_and_plot acs\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Modify-a-model","page":"API Documentation","title":"Modify a model","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"We list common transition attributes:","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"attribute interpretation\ntransPriority priority of a transition (influences resource allocation)\ntransProbOfSuccess probability that a transition terminates successfully\ntransCapacity maximum number of concurrent instances of the transition\ntransCycleTime duration of a transition's instance (adjusted by resource allocation)\ntransMaxLifeTime maximal duration of a transition's instance\ntransPostAction action to be executed once a transition's instance terminates\ntransName name of a transition","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"We list common species attributes:","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"attribute interpretation\nspecInitUncertainty uncertainty about variable's initial state (modelled as Gaussian standard deviation)\nspecInitVal initial value of a variable","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"Moreover, it is possible to specify the semantics of the \"rate\" term. By default, at each time step n ~ Poisson(rate * dt) instances of a given transition will be spawned. If you want to specify the rate in terms of a cycle time, you may want to use @ct(cycle_time), e.g., @ct(ex), A --> B, .... This is a shorthand for 1/ex, A --> B, ....","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"For deterministic \"rates\", use @per_step(ex). Here, ex evaluates to a deterministic number (ceiled to the nearest integer) of a transition's instances to spawn per a single integrator's step. However, note that in this case, the number doesn't scale with the step length! Moreover","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@add_species\n@aka\n@mode\n@name_transition","category":"page"},{"location":"index.html#ReactiveDynamics.@add_species","page":"API Documentation","title":"ReactiveDynamics.@add_species","text":"Add new species to a model.\n\nExamples\n\n@add_species acs S I R\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@aka","page":"API Documentation","title":"ReactiveDynamics.@aka","text":"Alias object name in an acs.\n\nDefault names\n\nname short name\nspecies S\ntransition T\naction A\nevent E\nparam P\nmeta M\n\nExamples\n\n@aka acs species=resource transition=reaction\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@mode","page":"API Documentation","title":"ReactiveDynamics.@mode","text":"Set species modality.\n\nSupported modalities\n\nnonblock\nconserved\nrate\n\nExamples\n\n@mode acs (r\"proj\\w+\", r\"experimental\\w+\") conserved\n@mode acs (S, I) conserved\n@mode acs S conserved\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@name_transition","page":"API Documentation","title":"ReactiveDynamics.@name_transition","text":"Set name of a transition in the model.\n\nExamples\n\n@name_transition acs 1=\"name\"\n@name_transition acs name=\"transition_name\"\n@name_transition acs \"name\"=\"transition_name\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Resource-costs","page":"API Documentation","title":"Resource costs","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@cost\n@valuation\n@reward","category":"page"},{"location":"index.html#ReactiveDynamics.@cost","page":"API Documentation","title":"ReactiveDynamics.@cost","text":"Set cost.\n\nExamples\n\n@cost model experimental1=2 experimental2=3\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@valuation","page":"API Documentation","title":"ReactiveDynamics.@valuation","text":"Set valuation.\n\nExamples\n\n@valuation model experimental1=2 experimental2=3\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@reward","page":"API Documentation","title":"ReactiveDynamics.@reward","text":"Set reward.\n\nExamples\n\n@reward model experimental1=2 experimental2=3\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Add-reactions","page":"API Documentation","title":"Add reactions","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@push\n@jump\n@periodic","category":"page"},{"location":"index.html#ReactiveDynamics.@push","page":"API Documentation","title":"ReactiveDynamics.@push","text":"Add reactions to an acset.\n\nExamples\n\n@push sir_acs β*S*I*tdecay(@time()) S+I --> 2I name=>SI2I\n@push sir_acs begin \n    ν*I, I --> R, name=>I2R\n    γ, R --> S, name=>R2S\nend\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@jump","page":"API Documentation","title":"ReactiveDynamics.@jump","text":"Add a jump process (with specified Poisson intensity per unit time step) to a model.\n\nExamples\n\n@jump acs λ Z += rand(Poisson(1.))\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@periodic","page":"API Documentation","title":"ReactiveDynamics.@periodic","text":"Add a periodic callback to a model.\n\nExamples\n\n@periodic acs 1. X += 1\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Set-initial-values,-uncertainty,-and-solver-arguments","page":"API Documentation","title":"Set initial values, uncertainty, and solver arguments","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@prob_init\n@prob_uncertainty\n@prob_params\n@prob_meta","category":"page"},{"location":"index.html#ReactiveDynamics.@prob_init","page":"API Documentation","title":"ReactiveDynamics.@prob_init","text":"Set initial values of species in an acset.\n\nExamples\n\n@prob_init acs X=1 Y=2 Z=h(α)\n@prob_init acs [1., 2., 3.]\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@prob_uncertainty","page":"API Documentation","title":"ReactiveDynamics.@prob_uncertainty","text":"Set uncertainty in initial values of species in an acset (stderr).\n\nExamples\n\n@prob_uncertainty acs X=.1 Y=.2\n@prob_uncertainty acs [.1, .2,]\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@prob_params","page":"API Documentation","title":"ReactiveDynamics.@prob_params","text":"Set parameter values in an acset.\n\nExamples\n\n@prob_params acs α=1. β=2.\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@prob_meta","page":"API Documentation","title":"ReactiveDynamics.@prob_meta","text":"Set model metadata (e.g. solver arguments)\n\nExamples\n\n@prob_meta acs tspan=(0, 100.) schedule=schedule_weighted!\n@prob_meta sir_acs tspan=250 tstep=1\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Model-unions","page":"API Documentation","title":"Model unions","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@join\n@equalize","category":"page"},{"location":"index.html#ReactiveDynamics.@join","page":"API Documentation","title":"ReactiveDynamics.@join","text":"@join models... [equalize...]\n\nPerforms join of models and identifies model variables, as specified.\n\nModel variables / parameter values and metadata are propagated; the last model takes precedence.\n\nExamples\n\n@join acs1 acs2 @catchall(A)=acs2.Z @catchall(XY) @catchall(B)\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@equalize","page":"API Documentation","title":"ReactiveDynamics.@equalize","text":"Identify (collapse) a set of species in a model.\n\nExamples\n\n@join acs acs1.A=acs2.A B=C\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Model-import-and-export","page":"API Documentation","title":"Model import and export","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@import_network\n@export_network","category":"page"},{"location":"index.html#ReactiveDynamics.@import_network","page":"API Documentation","title":"ReactiveDynamics.@import_network","text":"Import a model from a file: this can be either a single TOML file encoding the entire model, or a batch of CSV files (a root file and a number of files, each per a class of objects).\n\nSee tutorials/loadsave for an example.\n\nExamples\n\n@import_network \"model.toml\"\n@import_network \"csv/model.toml\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_network","page":"API Documentation","title":"ReactiveDynamics.@export_network","text":"Export model to a file: this can be either a single TOML file encoding the entire model, or a batch of CSV files (a root file and a number of files, each per a class of objects).\n\nSee tutorials/loadsave for an example.\n\nExamples\n\n@export_network acs \"acs_data.toml\" # as a TOML\n@export_network acs \"csv/model.csv\" # as a CSV\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Solution-import-and-export","page":"API Documentation","title":"Solution import and export","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@import_solution\n@export_solution_as_table\n@export_solution_as_csv\n@export_solution","category":"page"},{"location":"index.html#ReactiveDynamics.@import_solution","page":"API Documentation","title":"ReactiveDynamics.@import_solution","text":"@import_solution \"sol.jld2\"\n@import_solution \"sol.jld2\" sol\n\nImport a solution from a file.\n\nExamples\n\n@import_solution \"sir_acs_sol/serialized/sol.jld2\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_solution_as_table","page":"API Documentation","title":"ReactiveDynamics.@export_solution_as_table","text":"@export_solution_as_table sol\n\nExport a solution as a DataFrame.\n\nExamples\n\n@export_solution_as_table sol\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_solution_as_csv","page":"API Documentation","title":"ReactiveDynamics.@export_solution_as_csv","text":"@export_solution_as_csv sol\n@export_solution_as_csv sol \"sol.csv\"\n\nExport a solution to a file.\n\nExamples\n\n@export_solution_as_csv sol \"sol.csv\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_solution","page":"API Documentation","title":"ReactiveDynamics.@export_solution","text":"@export_solution sol\n@export_solution sol \"sol.jld2\"\n\nExport a solution to a file.\n\nExamples\n\n@export_solution sol \"sol.jdl2\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Problematize,sSolve,-and-plot","page":"API Documentation","title":"Problematize,sSolve, and plot","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@problematize\n@solve\n@plot","category":"page"},{"location":"index.html#ReactiveDynamics.@problematize","page":"API Documentation","title":"ReactiveDynamics.@problematize","text":"Convert a model to a DiscreteProblem. If passed a problem instance, return the instance.\n\nExamples\n\n@problematize acs tspan=1:100\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@solve","page":"API Documentation","title":"ReactiveDynamics.@solve","text":"Solve the problem. Solverargs passed at the calltime take precedence.\n\nExamples\n\n@solve prob\n@solve prob tspan=1:100\n@solve prob tspan=100 trajectories=20\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@plot","page":"API Documentation","title":"ReactiveDynamics.@plot","text":"Plot the solution (summary).\n\nExamples\n\n@plot sol plot_type=summary\n@plot sol plot_type=allocation # not supported for ensemble solutions!\n@plot sol plot_type=valuations # not supported for ensemble solutions!\n@plot sol plot_type=new_transitions # not supported for ensemble solutions!\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Optimization-and-fitting","page":"API Documentation","title":"Optimization and fitting","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@optimize\n@fit\n@fit_and_plot\n@build_solver","category":"page"},{"location":"index.html#ReactiveDynamics.@optimize","page":"API Documentation","title":"ReactiveDynamics.@optimize","text":"@optimize acset objective <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset and optimize given functional.\n\nObjective is an expression which may reference the model's variables and parameters, i.e., A+β. The values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model. The vector of free variables passed to the NLopt solver has the form [free_vars; free_params]; order of vars and params, respectively, is preserved. \n\nBy default, the functional is minimized. Specify objective=max to perform maximization. \n\nPropagates NLopt solver arguments; see NLopt documentation.\n\nExamples\n\n@optimize acs abs(A-B) A B=20. α=2. lower_bounds=0 upper_bounds=100\n@optimize acss abs(A-B) A B=20. α=2. upper_bounds=[200,300,400] maxeval=200 objective=min\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@fit","page":"API Documentation","title":"ReactiveDynamics.@fit","text":"@fit acset data_points time_steps empiric_variables <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset and fit initial values and parameters to empirical data.\n\nThe values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model. The vector of free variables passed to the NLopt solver has the form [free_vars; free_params]; order of vars and params, respectively, is preserved. \n\nPropagates NLopt solver arguments; see NLopt documentation.\n\nExamples\n\nt = [1, 50, 100]\ndata = [80 30 20]\n@fit acs data t vars=A B=20 A α # fit B, A, α; empirical data is for variable A\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@fit_and_plot","page":"API Documentation","title":"ReactiveDynamics.@fit_and_plot","text":"@fit acset data_points time_steps empiric_variables <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset, fit initial values and parameters to empirical data, and plot the result.\n\nThe values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model. The vector of free variables passed to the NLopt solver has the form [free_vars; free_params]; order of vars and params, respectively, is preserved. \n\nPropagates NLopt solver arguments; see NLopt documentation.\n\nExamples\n\nt = [1, 50, 100]\ndata = [80 30 20]\n@fit acs data t vars=A B=20 A α # fit B, A, α; empirical data is for variable A\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@build_solver","page":"API Documentation","title":"ReactiveDynamics.@build_solver","text":"@build_solver acset <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset and export a solution as a function of free vars and free parameters.\n\nExamples\n\nsolver = @build_solver acs S α β # function of variable S and parameters α, β\nsolver([S, α, β])\n\n\n\n\n\n","category":"macro"}]
+[{"location":"index.html#API-Documentation","page":"API Documentation","title":"API Documentation","text":"","category":"section"},{"location":"index.html#Create-a-model","page":"API Documentation","title":"Create a model","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@ReactionNetwork","category":"page"},{"location":"index.html#ReactiveDynamics.@ReactionNetwork","page":"API Documentation","title":"ReactiveDynamics.@ReactionNetwork","text":"Macro that takes an expression corresponding to a reaction network and outputs an instance of TheoryReactionNetwork that can be converted to a DiscreteProblem or solved directly.\n\nMost arrows accepted (both right, left, and bi-drectional arrows). Use 0 or ∅ for annihilation/creation to/from nothing.\n\nCustom functions and sampleable objects can be used as numeric parameters. Note that these have to be accessible from ReactiveDynamics's source code.\n\nExamples\n\nacs = @ReactionNetworkSchema begin\n    1.0, X ⟶ Y\n    1.0, X ⟶ Y, priority=>6., prob=>.7, capacity=>3.\n    1.0, ∅ --> (Poisson(.3γ)X, Poisson(.5)Y)\n    (XY > 100) && (XY -= 1)\nend\n@push acs 1.0 X ⟶ Y \n@prob_init acs X=1 Y=2 XY=α\n@prob_params acs γ=1 α=4\n@solve_and_plot acs\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Modify-a-model","page":"API Documentation","title":"Modify a model","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"We list common transition attributes:","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"attribute interpretation\ntransPriority priority of a transition (influences resource allocation)\ntransProbOfSuccess probability that a transition terminates successfully\ntransCapacity maximum number of concurrent instances of the transition\ntransCycleTime duration of a transition's instance (adjusted by resource allocation)\ntransMaxLifeTime maximal duration of a transition's instance\ntransPostAction action to be executed once a transition's instance terminates\ntransName name of a transition","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"We list common species attributes:","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"attribute interpretation\nspecInitUncertainty uncertainty about variable's initial state (modelled as Gaussian standard deviation)\nspecInitVal initial value of a variable","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"Moreover, it is possible to specify the semantics of the \"rate\" term. By default, at each time step n ~ Poisson(rate * dt) instances of a given transition will be spawned. If you want to specify the rate in terms of a cycle time, you may want to use @ct(cycle_time), e.g., @ct(ex), A --> B, .... This is a shorthand for 1/ex, A --> B, ....","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"For deterministic \"rates\", use @deterministic(ex). Here, ex evaluates to a deterministic number (ceiled to the nearest integer) of a transition's instances to spawn per a single integrator's step. However, note that in this case, the number doesn't scale with the step length! Moreover","category":"page"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@add_species\n@aka\n@mode\n@name_transition","category":"page"},{"location":"index.html#ReactiveDynamics.@add_species","page":"API Documentation","title":"ReactiveDynamics.@add_species","text":"Add new species to a model.\n\nExamples\n\n@add_species acs S I R\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@aka","page":"API Documentation","title":"ReactiveDynamics.@aka","text":"Alias object name in an acs.\n\nDefault names\n\nname short name\nspecies S\ntransition T\naction A\nevent E\nparam P\nmeta M\n\nExamples\n\n@aka acs species=resource transition=reaction\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@mode","page":"API Documentation","title":"ReactiveDynamics.@mode","text":"Set species modality.\n\nSupported modalities\n\nnonblock\nconserved\nrate\n\nExamples\n\n@mode acs (r\"proj\\w+\", r\"experimental\\w+\") conserved\n@mode acs (S, I) conserved\n@mode acs S conserved\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@name_transition","page":"API Documentation","title":"ReactiveDynamics.@name_transition","text":"Set name of a transition in the model.\n\nExamples\n\n@name_transition acs 1=\"name\"\n@name_transition acs name=\"transition_name\"\n@name_transition acs \"name\"=\"transition_name\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Resource-costs","page":"API Documentation","title":"Resource costs","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@cost\n@valuation\n@reward","category":"page"},{"location":"index.html#ReactiveDynamics.@cost","page":"API Documentation","title":"ReactiveDynamics.@cost","text":"Set cost.\n\nExamples\n\n@cost model experimental1=2 experimental2=3\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@valuation","page":"API Documentation","title":"ReactiveDynamics.@valuation","text":"Set valuation.\n\nExamples\n\n@valuation model experimental1=2 experimental2=3\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@reward","page":"API Documentation","title":"ReactiveDynamics.@reward","text":"Set reward.\n\nExamples\n\n@reward model experimental1=2 experimental2=3\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Add-reactions","page":"API Documentation","title":"Add reactions","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@push\n@jump\n@periodic","category":"page"},{"location":"index.html#ReactiveDynamics.@push","page":"API Documentation","title":"ReactiveDynamics.@push","text":"Add reactions to an acset.\n\nExamples\n\n@push sir_acs β*S*I*tdecay(@time()) S+I --> 2I name=>SI2I\n@push sir_acs begin \n    ν*I, I --> R, name=>I2R\n    γ, R --> S, name=>R2S\nend\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@jump","page":"API Documentation","title":"ReactiveDynamics.@jump","text":"Add a jump process (with specified Poisson intensity per unit time step) to a model.\n\nExamples\n\n@jump acs λ Z += rand(Poisson(1.))\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@periodic","page":"API Documentation","title":"ReactiveDynamics.@periodic","text":"Add a periodic callback to a model.\n\nExamples\n\n@periodic acs 1. X += 1\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Set-initial-values,-uncertainty,-and-solver-arguments","page":"API Documentation","title":"Set initial values, uncertainty, and solver arguments","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@prob_init\n@prob_uncertainty\n@prob_params\n@prob_meta","category":"page"},{"location":"index.html#ReactiveDynamics.@prob_init","page":"API Documentation","title":"ReactiveDynamics.@prob_init","text":"Set initial values of species in an acset.\n\nExamples\n\n@prob_init acs X=1 Y=2 Z=h(α)\n@prob_init acs [1., 2., 3.]\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@prob_uncertainty","page":"API Documentation","title":"ReactiveDynamics.@prob_uncertainty","text":"Set uncertainty in initial values of species in an acset (stderr).\n\nExamples\n\n@prob_uncertainty acs X=.1 Y=.2\n@prob_uncertainty acs [.1, .2,]\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@prob_params","page":"API Documentation","title":"ReactiveDynamics.@prob_params","text":"Set parameter values in an acset.\n\nExamples\n\n@prob_params acs α=1. β=2.\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@prob_meta","page":"API Documentation","title":"ReactiveDynamics.@prob_meta","text":"Set model metadata (e.g. solver arguments)\n\nExamples\n\n@prob_meta acs tspan=(0, 100.) schedule=schedule_weighted!\n@prob_meta sir_acs tspan=250 tstep=1\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Model-unions","page":"API Documentation","title":"Model unions","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@join\n@equalize","category":"page"},{"location":"index.html#ReactiveDynamics.@join","page":"API Documentation","title":"ReactiveDynamics.@join","text":"@join models... [equalize...]\n\nPerforms join of models and identifies model variables, as specified.\n\nModel variables / parameter values and metadata are propagated; the last model takes precedence.\n\nExamples\n\n@join acs1 acs2 @catchall(A)=acs2.Z @catchall(XY) @catchall(B)\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@equalize","page":"API Documentation","title":"ReactiveDynamics.@equalize","text":"Identify (collapse) a set of species in a model.\n\nExamples\n\n@join acs acs1.A=acs2.A B=C\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Model-import-and-export","page":"API Documentation","title":"Model import and export","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@import_network\n@export_network","category":"page"},{"location":"index.html#ReactiveDynamics.@import_network","page":"API Documentation","title":"ReactiveDynamics.@import_network","text":"Import a model from a file: this can be either a single TOML file encoding the entire model, or a batch of CSV files (a root file and a number of files, each per a class of objects).\n\nSee tutorials/loadsave for an example.\n\nExamples\n\n@import_network \"model.toml\"\n@import_network \"csv/model.toml\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_network","page":"API Documentation","title":"ReactiveDynamics.@export_network","text":"Export model to a file: this can be either a single TOML file encoding the entire model, or a batch of CSV files (a root file and a number of files, each per a class of objects).\n\nSee tutorials/loadsave for an example.\n\nExamples\n\n@export_network acs \"acs_data.toml\" # as a TOML\n@export_network acs \"csv/model.csv\" # as a CSV\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Solution-import-and-export","page":"API Documentation","title":"Solution import and export","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@import_solution\n@export_solution_as_table\n@export_solution_as_csv\n@export_solution","category":"page"},{"location":"index.html#ReactiveDynamics.@import_solution","page":"API Documentation","title":"ReactiveDynamics.@import_solution","text":"@import_solution \"sol.jld2\"\n@import_solution \"sol.jld2\" sol\n\nImport a solution from a file.\n\nExamples\n\n@import_solution \"sir_acs_sol/serialized/sol.jld2\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_solution_as_table","page":"API Documentation","title":"ReactiveDynamics.@export_solution_as_table","text":"@export_solution_as_table sol\n\nExport a solution as a DataFrame.\n\nExamples\n\n@export_solution_as_table sol\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_solution_as_csv","page":"API Documentation","title":"ReactiveDynamics.@export_solution_as_csv","text":"@export_solution_as_csv sol\n@export_solution_as_csv sol \"sol.csv\"\n\nExport a solution to a file.\n\nExamples\n\n@export_solution_as_csv sol \"sol.csv\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@export_solution","page":"API Documentation","title":"ReactiveDynamics.@export_solution","text":"@export_solution sol\n@export_solution sol \"sol.jld2\"\n\nExport a solution to a file.\n\nExamples\n\n@export_solution sol \"sol.jdl2\"\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Problematize,sSolve,-and-plot","page":"API Documentation","title":"Problematize,sSolve, and plot","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@problematize\n@solve\n@plot","category":"page"},{"location":"index.html#ReactiveDynamics.@problematize","page":"API Documentation","title":"ReactiveDynamics.@problematize","text":"Convert a model to a DiscreteProblem. If passed a problem instance, return the instance.\n\nExamples\n\n@problematize acs tspan=1:100\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@solve","page":"API Documentation","title":"ReactiveDynamics.@solve","text":"Solve the problem. Solverargs passed at the calltime take precedence.\n\nExamples\n\n@solve prob\n@solve prob tspan=1:100\n@solve prob tspan=100 trajectories=20\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@plot","page":"API Documentation","title":"ReactiveDynamics.@plot","text":"Plot the solution (summary).\n\nExamples\n\n@plot sol plot_type=summary\n@plot sol plot_type=allocation # not supported for ensemble solutions!\n@plot sol plot_type=valuations # not supported for ensemble solutions!\n@plot sol plot_type=new_transitions # not supported for ensemble solutions!\n\n\n\n\n\n","category":"macro"},{"location":"index.html#Optimization-and-fitting","page":"API Documentation","title":"Optimization and fitting","text":"","category":"section"},{"location":"index.html","page":"API Documentation","title":"API Documentation","text":"@optimize\n@fit\n@fit_and_plot\n@build_solver","category":"page"},{"location":"index.html#ReactiveDynamics.@optimize","page":"API Documentation","title":"ReactiveDynamics.@optimize","text":"@optimize acset objective <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset and optimize given functional.\n\nObjective is an expression which may reference the model's variables and parameters, i.e., A+β. The values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model. The vector of free variables passed to the NLopt solver has the form [free_vars; free_params]; order of vars and params, respectively, is preserved. \n\nBy default, the functional is minimized. Specify objective=max to perform maximization. \n\nPropagates NLopt solver arguments; see NLopt documentation.\n\nExamples\n\n@optimize acs abs(A-B) A B=20. α=2. lower_bounds=0 upper_bounds=100\n@optimize acss abs(A-B) A B=20. α=2. upper_bounds=[200,300,400] maxeval=200 objective=min\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@fit","page":"API Documentation","title":"ReactiveDynamics.@fit","text":"@fit acset data_points time_steps empiric_variables <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset and fit initial values and parameters to empirical data.\n\nThe values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model. The vector of free variables passed to the NLopt solver has the form [free_vars; free_params]; order of vars and params, respectively, is preserved. \n\nPropagates NLopt solver arguments; see NLopt documentation.\n\nExamples\n\nt = [1, 50, 100]\ndata = [80 30 20]\n@fit acs data t vars=A B=20 A α # fit B, A, α; empirical data is for variable A\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@fit_and_plot","page":"API Documentation","title":"ReactiveDynamics.@fit_and_plot","text":"@fit acset data_points time_steps empiric_variables <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset, fit initial values and parameters to empirical data, and plot the result.\n\nThe values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model. The vector of free variables passed to the NLopt solver has the form [free_vars; free_params]; order of vars and params, respectively, is preserved. \n\nPropagates NLopt solver arguments; see NLopt documentation.\n\nExamples\n\nt = [1, 50, 100]\ndata = [80 30 20]\n@fit acs data t vars=A B=20 A α # fit B, A, α; empirical data is for variable A\n\n\n\n\n\n","category":"macro"},{"location":"index.html#ReactiveDynamics.@build_solver","page":"API Documentation","title":"ReactiveDynamics.@build_solver","text":"@build_solver acset <free_var=[init_val]>... <free_prm=[init_val]>... opts...\n\nTake an acset and export a solution as a function of free vars and free parameters.\n\nExamples\n\nsolver = @build_solver acs S α β # function of variable S and parameters α, β\nsolver([S, α, β])\n\n\n\n\n\n","category":"macro"}]
 }
diff --git a/docs/src/index.md b/docs/src/index.md
index 87ec4e7..7104f87 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -28,7 +28,7 @@ We list common species attributes:
 
 Moreover, it is possible to specify the semantics of the "rate" term. By default, at each time step `n ~ Poisson(rate * dt)` instances of a given transition will be spawned. If you want to specify the rate in terms of a cycle time, you may want to use `@ct(cycle_time)`, e.g., `@ct(ex), A --> B, ...`. This is a shorthand for `1/ex, A --> B, ...`.
 
-For deterministic "rates", use `@per_step(ex)`. Here, `ex` evaluates to a deterministic number (ceiled to the nearest integer) of a transition's instances to spawn per a single integrator's step. However, note that in this case, the number doesn't scale with the step length! Moreover
+For deterministic "rates", use `@deterministic(ex)`. Here, `ex` evaluates to a deterministic number (ceiled to the nearest integer) of a transition's instances to spawn per a single integrator's step. However, note that in this case, the number doesn't scale with the step length! Moreover
 
 ```@docs
 @add_species
diff --git a/src/ReactiveDynamics.jl b/src/ReactiveDynamics.jl
index c10b7a0..7d4991d 100644
--- a/src/ReactiveDynamics.jl
+++ b/src/ReactiveDynamics.jl
@@ -104,7 +104,7 @@ Base.convert(::Type{FoldedObservable}, ex::String) = eval(Meta.parse(ex))
 prettynames = Dict(
     :transRate => [:rate],
     :specInitUncertainty => [:uncertainty, :stoch, :stochasticity],
-    :transPreAction => [:preAction, :pre],
+    :transPreAction => [:preAction, :action, :pre],
     :transPostAction => [:postAction, :post],
     :transName => [:name, :interpretation],
     :transPriority => [:priority],
diff --git a/src/compilers.jl b/src/compilers.jl
index eb651ba..49ecffd 100644
--- a/src/compilers.jl
+++ b/src/compilers.jl
@@ -64,8 +64,7 @@ function recursively_expand_dots_in_ex!(ex, vars)
     return ex
 end
 
-reserved_names =
-    [:t, :state, :obs, :resample, :solverarg, :take, :log, :periodic, :set_params]
+reserved_names = [:t, :obs, :resample, :solverarg, :take, :log, :periodic, :set_params]
 
 function escape_ref(ex, species)
     return if ex isa Symbol
@@ -92,20 +91,16 @@ function wrap_expr(fex, species_names, prm_names, varmap)
         let
         end
     )
+
     # expression walking (MacroTools): visit each expression, subsitute with the body's return value
     fex = prewalk(fex) do x
         # here we convert the query metalanguage: @t() -> time(state) etc. 
         if isexpr(x, :macrocall) && (macroname(x) ∈ reserved_names)
             Expr(:call, macroname(x), :state, x.args[3:end]...)
-        else
-            x
-        end
-    end
-
-    fex = prewalk(fex) do x
-        # here we convert the query metalanguage: @t() -> time(state) etc. 
-        if isexpr(x, :macrocall) && (macroname(x) == :transition)
+        elseif isexpr(x, :macrocall) && (macroname(x) == :transition)
             :transition
+        elseif isexpr(x, :macrocall) && (macroname(x) == :state)
+            :state
         else
             x
         end
@@ -148,7 +143,7 @@ function skip_compile(attr)
 end
 
 function compile_attrs(acs::ReactionNetworkSchema, structured_token)
-    species_names = setdiff(collect(acs[:, :specName]), structured_token)
+    species_names = collect(acs[:, :specName])#setdiff(collect(acs[:, :specName]), structured_token)
 
     prm_names = collect(acs[:, :prmName])
     varmap = Dict([name => :(state.u[$i]) for (i, name) in enumerate(species_names)])
diff --git a/src/interface/create.jl b/src/interface/create.jl
index 099ca1e..faf5bdf 100644
--- a/src/interface/create.jl
+++ b/src/interface/create.jl
@@ -146,7 +146,7 @@ function recursively_expand_actions!(evs, condex, event)
 end
 
 function expand_rate(rate)
-    rate = if !(isexpr(rate, :macrocall) && (macroname(rate) == :per_step))
+    rate = if !(isexpr(rate, :macrocall) && (macroname(rate) == :deterministic))
         :(rand(Poisson(max(state.dt * $rate, 0))))
     else
         rate.args[3]
@@ -271,7 +271,7 @@ function prune_reaction_line!(pcs, reactants, line)
     return line
 end
 
-function recursively_find_reactants!(reactants, pcs, ex::SampleableValues)
+function recursively_find_reactants!(reactants, pcs, ex)
     if typeof(ex) != Expr || isexpr(ex, :.) || (ex.head == :escape)
         if (ex == 0 || in(ex, empty_set))
             return :∅
@@ -293,8 +293,11 @@ function recursively_find_reactants!(reactants, pcs, ex::SampleableValues)
                 isexpr(ex.args[i], :tuple) ? ex.args[i].args[2] : ex.args[i],
             )
         end
+    elseif isexpr(ex, :macrocall) && macroname(ex) ∈ [:structured, :move]
+        return ex
     elseif isexpr(ex, :macrocall)
-        recursively_find_reactants!(reactants, pcs, ex.args[3])
+        pass_value = ex.args[3] isa QuoteNode ? ex.args[3].value : ex.args[3]
+        recursively_find_reactants!(reactants, pcs, pass_value)
     elseif isexpr(ex, :call)
         push!(reactants, ex.args[1])
     else
diff --git a/src/solvers.jl b/src/solvers.jl
index 9bf982c..46d3e6a 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -332,19 +332,31 @@ end
 
 function structured_rhs(expr::Expr, state, transition)
     if isexpr(expr, :macrocall) && macroname(expr) == :structured
-        expr = quote
-            token = $(expr.args[end-1])
-            species = $(expr.args[end])
+        if length(expr.args) == 3
+            expr = quote
+                return $(expr.args[end])
+            end
+            # write docs
+            token = context_eval(state, transition, state.wrap_fun(expr))
 
-            return token, species
-        end
+            entangle!(getagent(state, "structured"), token)
 
-        token, species = context_eval(state, transition, state.wrap_fun(expr))
-        set_species!(token, Symbol(species))
+            return token, get_species(token)
+        else
+            expr = quote
+                token = $(expr.args[end-1])
+                species = $(expr.args[end])
 
-        entangle!(getagent(state, "structured"), token)
+                return token, species
+            end
+            # write docs
+            token, species = context_eval(state, transition, state.wrap_fun(expr))
+            set_species!(token, Symbol(species))
 
-        return token, get_species(token)
+            entangle!(getagent(state, "structured"), token)
+
+            return token, get_species(token)
+        end
     elseif isexpr(expr, :macrocall) && macroname(expr) == :move
         expr = quote
             species_from = $(expr.args[end-1])
@@ -528,6 +540,8 @@ function ReactionNetworkProblem(
     structured_token_names =
         acs[filter(i -> acs[i, :specStructured], 1:nparts(acs, :S)), :specName]
 
+    println(acs[:, :specName])
+    println(structured_token_names)
     attrs, transitions, wrap_fun = compile_attrs(acs, structured_token_names)
     transition_recipes = transitions
     u0_init = zeros(nparts(acs, :S))
@@ -583,7 +597,7 @@ function ReactionNetworkProblem(
 
     entangle!(network, FreeAgent("structured"))
 
-    save!(network)
+    # save!(network)
 
     return network
 end
@@ -612,6 +626,11 @@ function update_u_structured!(state)
 end
 
 function AlgebraicAgents._step!(state::ReactionNetworkProblem)
+    update_u_structured!(state)
+    if isempty(state.sol)
+        save!(state)
+    end
+
     free_blocked_species!(state)
     update_u_structured!(state)
     update_observables(state)
@@ -632,8 +651,11 @@ function AlgebraicAgents._step!(state::ReactionNetworkProblem)
         ),
     )
 
+    state.t += state.dt
+
     save!(state)
-    return state.t += state.dt
+
+    return state.t
 end
 
 function AlgebraicAgents._projected_to(state::ReactionNetworkProblem)
diff --git a/tutorial/agents-integration/Project.toml b/tutorial/agents-integration/Project.toml
index ef794fc..dc9c60a 100644
--- a/tutorial/agents-integration/Project.toml
+++ b/tutorial/agents-integration/Project.toml
@@ -1,9 +1,16 @@
 [deps]
 AlgebraicAgents = "f6eb0ae3-10fa-40e6-88dd-9006ba45093a"
+CEEDesigns = "e939450b-799e-4198-a5f5-3f2f7fb1c671"
+Copulas = "ae264745-0b69-425e-9d9d-cf662c5eec93"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+MCTS = "e12ccd36-dcad-5f33-8774-9175229e7b33"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+POMDPTools = "7588e00f-9cae-40de-98dc-e0c70c48cdd7"
+POMDPs = "a93abf59-7444-517b-a68a-c42f96afdd7d"
 PlotGraphviz = "78a92bc3-407c-4e2f-aae5-75bb47a6fe36"
 ReactiveDynamics = "c7456e7d-545a-4b79-91ea-6e93d96dd4d4"
+ScientificTypes = "321657f4-b219-11e9-178b-2701a2544e81"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
diff --git a/tutorial/agents-integration/agents.jl b/tutorial/agents-integration/agents.jl
index b021720..6b6bc4d 100644
--- a/tutorial/agents-integration/agents.jl
+++ b/tutorial/agents-integration/agents.jl
@@ -24,7 +24,7 @@ network = @ReactionNetworkSchema
     ρ2, ∅ --> R2
 
     # Generate "Molecule 1" (where the integer corresponds to a "state" of, e.g., experimental triage).
-    ρ3, ∅ --> M1(@t(), rand(4))
+    ρ3, ∅ --> @structured(M1(@t(), rand(4)))
 
     # Based on properties of particular "structured agent" assigned to the transition,
     # we can update the attributes of the instance of a transition (such as probability of success).
@@ -33,28 +33,26 @@ network = @ReactionNetworkSchema
     # Update transition probability based on properties of "M1," 
     # which was assigned as a "resource" to the transition.
     ρ4,
-    R1 + M1 --> M2(@t(), rand(4)),
+    R1 + SM1 --> @structured(M2(@t(), rand(4))),
     preAction => update_prob_transition(state, transition)
 
-    # R2 + M1 --> M2(@t(), rand(4))
-    5.0, R2 + M1 --> @structured(M2(@t(), rand(4)), :A)
+    5.0, R2 + SM1 --> @structured(M2(@t(), rand(4)), :A)
     1.0, R2 + A --> @structured(M2(@t(), rand(4)), f_species(@transition))
 
-    1.0, R2 + M1 --> @move(M1, :M2)
-    1.0, R2 + M1 --> @move(M1, :C)
+    2.0, R2 + SM1 --> @move(:SM1, :C)
 end
 
 @prob_init network R1 = 10 R2 = 15
 
 # As for structured agents, we will need to instantiate the instances
 # and add them to the instance of a network. But first, we still need to define these types.
-@prob_init network M1 = 2 M2 = 0
+@prob_init network SM1 = 2 SM2 = 0
 
 @prob_params network ρ1 = 2 ρ2 = 1 ρ3 = 3 ρ4 = 4
 
 @prob_meta network tspan = 100 dt = 1.0
 
-# We use `@structured` macro, which is a convenience wrapper around `@aagent`),
+# We use `@structured` macro, which is a convenience wrapper around `@aagent`,
 # defined in ReactiveDynamics.jl
 @structured_token network struct M1
     descriptor::Any
@@ -68,7 +66,7 @@ using Random
 using ReactiveDynamics: M1
 
 function ReactiveDynamics.M1(time, descriptor)
-    return M1("M1" * randstring(4), :M1, nothing, [], descriptor, time)
+    return M1("M1" * randstring(4), :SM1, nothing, [], descriptor, time)
 end
 
 # We define the function which updates the transition probability.
@@ -94,12 +92,12 @@ end
         time_created::Any
     end
 
-    using Random
-    M2(time, descriptor) = M2("M2" * randstring(4), :M2, nothing, [], descriptor, time)
+    using Random: randstring
+    M2(time, descriptor) = M2("M2" * randstring(4), :SM2, nothing, [], descriptor, time)
 end
 
 # Let the network know that the species is structured.
-for species in [:M1, :M2, :A, :B, :C]
+for species in [:SM1, :SM2, :A, :B, :C]
     ReactiveDynamics.register_structured_species!(network, species)
 end
 
diff --git a/tutorial/agents-integration/ceed_integration.jl b/tutorial/agents-integration/ceed_integration.jl
new file mode 100644
index 0000000..18a29fe
--- /dev/null
+++ b/tutorial/agents-integration/ceed_integration.jl
@@ -0,0 +1,236 @@
+using CEEDesigns, CEEDesigns.GenerativeDesigns
+using DataFrames
+using ScientificTypes
+using Statistics, Copulas
+import POMDPs, POMDPTools, MCTS
+
+import Distributions
+
+# ----- Experimental Setup -----
+
+# We generate a synthetic dataset.
+# This is taken from https://github.com/Merck/CEEDesigns.jl/blob/34588ae0e5563cb93f6818e3a9c8b3a77c5e3c47/tutorials/SimpleGenerative.jl
+
+include("experimental_setup.jl")
+
+# ----- Get a sampling function -----
+
+(; sampler, uncertainty, weights) = DistanceBased(
+    data;
+    target = "y",
+    uncertainty = Variance(),
+    similarity = GenerativeDesigns.Exponential(; λ = 5),
+);
+
+# ----- Set up a reaction network -----
+
+#=
+Pkg.activate(".")
+Pkg.develop(path = "../..")
+=#
+
+using ReactiveDynamics
+
+# Set up parameters that will be used to define a network.
+
+# Experiments and costs
+features_experiments = Dict(["x$i" => "e$i" for i = 1:4])
+
+experiments_costs = Dict([
+    features_experiments[e] => (i, i) => [e] for (i, e) in enumerate(names(data)[1:4])
+])
+
+experiments_costs["ey"] = (100, 100) => ["y"]
+
+# Experimental resources
+experimental_resources = [:ER1, :ER2, :ER3]
+resources_quantities = [rand(1:3, length(experimental_resources)) for _ = 1:5]
+
+# "Compound," which is a structured token.
+# We use `@structured` macro, which is a convenience wrapper around `@aagent`,
+# defined in ReactiveDynamics.jl
+@register begin
+    @aagent BaseStructuredToken AbstractStructuredToken struct Compound
+        state::Any
+        history::Vector{Symbol}
+    end
+
+    get_cmpds = function (transition::Transition) end
+    run_experiment = function (agent::Compound, experiment::Symbol, rng) end
+    assign_to_places = function (state::ReactionNetworkProblem, threshold) end
+end
+
+# Provide a constructor for `Compound` and define functions that will
+# "execute" the experiments.
+import ReactiveDynamics: Compound, get_cmpds, assign_to_places, sample
+using ReactiveDynamics: ReactionNetworkProblem, Transition
+
+function Compound(id::AbstractString, predictions::Dict)
+    state = State((Evidence(predictions...), Tuple(zeros(2))))
+
+    return Compound("Compound $id", :pool, nothing, [], state, String[])
+end
+
+using Random: default_rng
+
+ReactiveDynamics.get_cmpds = function (transition::Transition)
+    if !isnothing(transition) && !isempty(transition.bound_structured_agents)
+        agent_ix = findall(x -> x isa Compound, transition.bound_structured_agents)
+
+        return transition.bound_structured_agents[agent_ix]
+    end
+end
+
+ReactiveDynamics.run_experiment =
+    function (agents::Vector, experiment::Symbol, rng = default_rng())
+        println("running experiment $experiment")
+        for agent in agents
+            push!(agent.history, experiment)
+            experiment = String(experiment)
+
+            observation = sampler(
+                agent.state.evidence,
+                getindex(experiments_costs[experiment], 2),
+                rng,
+            )
+
+            agent.state =
+                merge(agent.state, observation, first(experiments_costs[experiment]))
+        end
+        return agents
+    end
+
+transitions_experiments = String[]
+for i = 1:5
+    experiment = i < 5 ? "e$i" : "ey"
+
+    resource_part = join(
+        [
+            "$(resources_quantities[i][res_i]) * $res" for
+            (res_i, res) in enumerate(experimental_resources)
+        ],
+        " + ",
+    )
+
+    push!(
+        transitions_experiments,
+        """@deterministic($experiment), $experiment + $resource_part --> @move(:$experiment, :pool),
+            action => run_experiment(get_cmpds(@transition), :$experiment)
+        """,
+    )
+end
+
+# Set up a reaction network
+
+network = @ReactionNetworkSchema
+
+# Resource generation part
+@push network begin
+    p1, ∅ --> ER1
+    p2, ∅ --> ER2
+    p3, ∅ --> ER3
+end
+
+@prob_init network ER1 = 1200 ER2 = 1500 ER3 = 1300
+
+@prob_params network p1 = 1 p2 = 1 p3 = 1
+
+# Experiments part
+for species in union(Symbol.(keys(experiments_costs)), [:pool])
+    ReactiveDynamics.register_structured_species!(network, species)
+end
+
+str_network_def = """
+    begin
+        @push network begin
+            $(join(transitions_experiments, '\n'))
+        end
+    end 
+"""
+
+using MacroTools: striplines
+expr_network_def = striplines(Meta.parseall(str_network_def))
+
+eval(expr_network_def)
+
+@prob_meta network tspan = 100
+
+problem = ReactionNetworkProblem(network)
+
+using Random: randstring
+
+# Simplified setup: we assume that compounds are already assigned
+# to the "experimental" places.
+for _ = 1:10
+    cmpd = Compound(randstring(4), Dict())
+    cmpd.species = Symbol("e$(rand(1:4))")
+
+    add_structured_token!(problem, cmpd)
+end
+
+simulate(problem, 10)
+
+# To allow for dynamic assignment of compounds to places, we need to create an agent
+# that will move the agent to the corresponding places.
+# This can be expressed in two ways:
+# - Create an "algebraic agent" which will modify the reaction network's state,
+# - Create a "placeholder" transition which will run the mutating function.
+
+# In either case, we need to define the function that will facilitate the assignments.
+evidence = Evidence()
+
+solver = GenerativeDesigns.DPWSolver(; n_iterations = 500, tree_in_info = true)
+repetitions = 5
+mdp_options = (; max_parallel = 1, discount = 1.0, costs_tradeoff = (0.5, 0.5))
+
+ReactiveDynamics.assign_to_places =
+    function (state::ReactionNetworkProblem, threshold = 0.1)
+        compounds = filter(
+            x -> ReactiveDynamics.get_species(x) == :pool,
+            collect(values(inners(inners(state)["structured"]))),
+        )
+
+        for cmpd in compounds
+            e = get_next_experiment(cmpd.state.evidence, threshold)
+            if !isnothing(e)
+                cmpd.species = first(e)
+            end
+        end
+    end
+
+function get_next_experiment(evidence::Evidence, threshold = 0.1)
+    design = efficient_design(
+        experiments_costs;
+        sampler = sampler,
+        uncertainty = uncertainty,
+        threshold,
+        evidence = evidence,
+        solver = solver,
+        repetitions = repetitions,
+        mdp_options = mdp_options,
+    )
+
+    if !haskey(design[2], :arrangement) || isempty(design[2].arrangement)
+        return nothing
+    else
+        return first(design[2].arrangement)
+    end
+end
+
+@push network begin
+    @deterministic(2), ∅ --> @structured(Compound(randstring(4), Dict()))
+    @deterministic(1), ∅ --> ER1, preAction => assign_to_places(@state)
+end
+
+problem = ReactionNetworkProblem(network)
+
+# Simplified setup: we assume that compounds are already assigned
+# to the "experimental" places.
+for _ = 1:10
+    cmpd = Compound(randstring(4), Dict())
+    cmpd.species = Symbol("e$(rand(1:4))")
+
+    add_structured_token!(problem, cmpd)
+end
+
+simulate(problem, 10)
diff --git a/tutorial/agents-integration/experimental_setup.jl b/tutorial/agents-integration/experimental_setup.jl
new file mode 100644
index 0000000..cb5c4d0
--- /dev/null
+++ b/tutorial/agents-integration/experimental_setup.jl
@@ -0,0 +1,56 @@
+# This is taken from https://github.com/Merck/CEEDesigns.jl/blob/34588ae0e5563cb93f6818e3a9c8b3a77c5e3c47/tutorials/SimpleGenerative.jl
+
+make_friedman3 = function (U, noise = 0, friedman3 = true)
+    size(U, 2) == 4 || error("input U must have 4 columns, has $(size(U,2))")
+    n = size(U, 1)
+    X = DataFrame(zeros(Float64, n, 4), :auto)
+    for i = 1:4
+        X[:, i] .= U[:, i]
+    end
+    ϵ = noise > 0 ? rand(Distributions.Normal(0, noise), size(X, 1)) : 0
+    if friedman3
+        X.y = @. atan((X[:, 2] * X[:, 3] - 1 / (X[:, 2] * X[:, 4])) / X[:, 1]) + ϵ
+    else
+        ## Friedman #2
+        X.y = @. (X[:, 1]^2 + (X[:, 2] * X[:, 3] - 1 / (X[:, 2] * X[:, 4]))^2)^0.5 + ϵ
+    end
+    return X
+end
+
+p12, p13, p14, p23, p24, p34 = 0.8, 0.5, 0.3, 0.5, 0.25, 0.4
+Σ = [
+    1 p12 p13 p14
+    p12 1 p23 p24
+    p13 p23 1 p34
+    p14 p24 p34 1
+]
+
+X1 = Distributions.Uniform(0, 100)
+X2 = Distributions.Uniform(40 * π, 560 * π)
+X3 = Distributions.Uniform(0, 1)
+X4 = Distributions.Uniform(1, 11)
+
+C = GaussianCopula(Σ)
+D = SklarDist(C, (X1, X2, X3, X4))
+
+X = rand(D, 1000)
+
+data = make_friedman3(transpose(X), 0.01)
+
+data[1:10, :]
+
+# We can check that the empirical correlation is roughly the same as the specified theoretical values: 
+
+cor(Matrix(data[:, Not(:y)]))
+
+# We ensure that our algorithms know that we have provided data of specified types. 
+
+types = Dict(
+    :x1 => ScientificTypes.Continuous,
+    :x2 => ScientificTypes.Continuous,
+    :x3 => ScientificTypes.Continuous,
+    :x4 => ScientificTypes.Continuous,
+    :y => ScientificTypes.Continuous,
+)
+
+data = coerce(data, types);
diff --git a/tutorial/toy_pharma_model.jl b/tutorial/toy_pharma_model.jl
index e536d19..c308fe7 100644
--- a/tutorial/toy_pharma_model.jl
+++ b/tutorial/toy_pharma_model.jl
@@ -45,19 +45,19 @@ draw(sol)
 
 # model dynamics
 toy_pharma_model = @ReactionNetworkSchema begin
-    @per_step(α(candidate_compound, marketed_drug, κ)),
+    @deterministic(α(candidate_compound, marketed_drug, κ)),
     3 * @conserved(scientist) + @rate(budget) --> candidate_compound,
     name => discovery,
     probability => 0.3,
     cycletime => 10.0,
     priority => 0.5
-    @per_step(β(candidate_compound, marketed_drug)),
+    @deterministic(β(candidate_compound, marketed_drug)),
     candidate_compound + 5 * @conserved(scientist) + 2 * @rate(budget) -->
     marketed_drug + 5 * budget,
     name => dx2market,
     probability => 0.5 + 0.001 * @t(),
     cycletime => 4
-    @per_step(γ * marketed_drug), marketed_drug --> ∅, name => drug_killed
+    @deterministic(γ * marketed_drug), marketed_drug --> ∅, name => drug_killed
 end
 
 @periodic toy_pharma_model 0.0 budget += 11 * marketed_drug