diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
new file mode 100644
index 0000000..de7ac75
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -0,0 +1,27 @@
+---
+name: Bug report
+about: Create a report to help us improve
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+**Describe the bug**
+A clear and concise description of what the bug is.
+
+**Minimal Working Example**
+Please provide a piece of code that leads to the bug you encounter.
+
+If the code is **runnable**, it will help us identify the problem faster.
+
+**Package versions**
+
+Please provide the versions you use. To do this, run the code:
+```julia
+using Pkg
+Pkg.status([
+    "Package1", "Package2"]; # etc.
+    mode = PKGMODE_MANIFEST
+)
+```
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
new file mode 100644
index 0000000..91bb1de
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/feature_request.md
@@ -0,0 +1,14 @@
+---
+name: Feature request
+about: Suggest an idea for this project
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+**Describe the feature you'd like to have**
+
+**Cite scientific papers related to the feature/algorithm**
+
+**If possible, sketch out an implementation strategy**
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
new file mode 100644
index 0000000..6c2390a
--- /dev/null
+++ b/.github/workflows/CompatHelper.yml
@@ -0,0 +1,24 @@
+name: CompatHelper
+
+on:
+  schedule:
+    - cron: '00 * * * *'
+
+jobs:
+  CompatHelper:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1]
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
new file mode 100644
index 0000000..623860f
--- /dev/null
+++ b/.github/workflows/TagBot.yml
@@ -0,0 +1,15 @@
+name: TagBot
+on:
+  issue_comment:
+    types:
+      - created
+  workflow_dispatch:
+jobs:
+  TagBot:
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          ssh: ${{ secrets.DOCUMENTER_KEY }}
\ No newline at end of file
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..6278b77
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,53 @@
+name: CI
+on:
+  pull_request:
+    branches:
+      - main
+      - '**' # matches every branch
+  push:
+    branches:
+      - main
+    tags: '*'
+jobs:
+  test:
+    name: Tests, Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1'
+        os:  [ubuntu-latest] # adjust according to need, e.g. os: [ubuntu-latest] if testing only on linux
+        arch:
+          - x64
+    steps:
+      # Cancel ongoing CI test runs if pushing to branch again before the previous tests
+      # have finished
+      - name: Cancel ongoing test runs for previous commits
+        uses: styfle/cancel-workflow-action@0.6.0
+        with:
+          access_token: ${{ github.token }}
+
+      # Do tests
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
diff --git a/.github/workflows/doccleanup.yml b/.github/workflows/doccleanup.yml
new file mode 100644
index 0000000..22d9043
--- /dev/null
+++ b/.github/workflows/doccleanup.yml
@@ -0,0 +1,26 @@
+name: Doc Preview Cleanup
+
+on:
+  pull_request:
+    types: [closed]
+
+jobs:
+  doc-preview-cleanup:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout gh-pages branch
+        uses: actions/checkout@v2
+        with:
+          ref: gh-pages
+      - name: Delete preview and history + push changes
+        run: |
+            if [ -d "previews/PR$PRNUM" ]; then
+              git config user.name "Documenter.jl"
+              git config user.email "documenter@juliadocs.github.io"
+              git rm -rf "previews/PR$PRNUM"
+              git commit -m "delete preview"
+              git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree})
+              git push --force origin gh-pages-new:gh-pages
+            fi
+        env:
+            PRNUM: ${{ github.event.number }}
\ No newline at end of file
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
new file mode 100644
index 0000000..33e080c
--- /dev/null
+++ b/.github/workflows/documentation.yml
@@ -0,0 +1,26 @@
+name: Documentation
+
+on:
+  push:
+    branches:
+      - main
+    tags: '*'
+  pull_request:
+
+jobs:
+  build:
+    permissions:
+      contents: write
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: '1'
+      - name: Install dependencies
+        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
+        run: julia --project=docs/ docs/make.jl
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..32639ac
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,12 @@
+name = "ProcessBasedModelling"
+uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb"
+authors = ["Datseris <datseris.george@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+
+[compat]
+ModelingToolkit = "8.73"
+Reexport = "1.2"
\ No newline at end of file
diff --git a/README.md b/README.md
index d52f7b2..8407f13 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,38 @@
-# ProcessBasedModelling.jl
\ No newline at end of file
+# ProcessBasedModelling.jl
+
+[![docsdev](https://img.shields.io/badge/docs-dev-lightblue.svg)](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/ProcessBasedModelling/dev/)
+[![docsstable](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/ProcessBasedModelling/stable/)
+[![Paper](https://img.shields.io/badge/Cite-DOI:10.1063/5.0159675-purple)](https://arxiv.org/abs/2304.12786)
+[![CI](https://github.com/JuliaDynamics/ProcessBasedModelling.jl/workflows/CI/badge.svg)](https://github.com/JuliaDynamics/ProcessBasedModelling.jl/actions?query=workflow%3ACI)
+[![codecov](https://codecov.io/gh/JuliaDynamics/ProcessBasedModelling.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDynamics/ProcessBasedModelling.jl)
+[![Package Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/ProcessBasedModelling)](https://pkgs.genieframework.com?packages=ProcessBasedModelling)
+
+ProcessBasedModelling.jl is an extension to [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) (MTK).
+It is an alternative framework to MTK's
+native component-based modelling, but, instead of components, there are "processes".
+This modelling approach is useful in the modelling of physical/biological/whatever systems, where each variable corresponds to particular physical concept or observable and there are few (or any) duplicate variables to make the definition of MTK "factories" worthwhile.
+In many modelling scenarios this follows the line of reasoning in the researcher's head.
+
+Beyond this reasoning style, the biggest strength of ProcessBasedModelling.jl is the informative errors it provides regarding incorrect/incomplete equations. When building the MTK model via ProcessBasedModelling.jl the user provides a vector of "processes": equations or custom types that have a well defined and single left-hand-side variable.
+This allows ProcessBasedModelling.jl to achieve the following:
+
+1. Iterates over the processes and collects _new_ variables that have been introduced by a provided process but do not themselves have a process assigned to them.
+2. For these collected "process-less" variables:
+   - If there is a default process defined, it incorporates this one into the model
+   - If there is no default process but the variable has a default value, it equates the variable to a _parameter_ that has the same default value and throws an informative warning.
+   - Else, it throws an informative error saying exactly which originally provided variable introduced this new "process-less" variable.
+3. Also throws an informative error if a variable has two processes assigned to it (by mistake).
+
+In our experience, and as we also highlight explicitly in the online documentation, this approach typically yields simpler, less ambiguous and more targeted warning/error messages than the native MTK one's, leading to faster identification and resolution of the problems with the composed equations.
+
+ProcessBasedModelling.jl is particularly suited for developing a model about a physical/biological/whatever system and being able to try various physical "rules" (couplings, feedbacks, mechanisms, ...) for a given physical observable efficiently.
+This means switching arbitrarily between different processes that correspond to the same variable.
+Hence, the target application of ProcessBasedModelling.jl is to be a framework to develop field-specific libraries that offer predefined processes without themselves relying on the existence of context-specific predefined components. An example usage is in [EnergyBalanceModels.jl](https://github.com/JuliaDynamics/EnergyBalanceModels.jl).
+
+Besides the informative errors, ProcessBasedModelling.jl also
+
+1. Provides a couple of common process subtypes out of the box to accelerate development of field-specific libraries.
+2. Makes named MTK variables and parameters automatically, corresponding to parameters introduced by the by-default provided processes. This typically leads to intuitive names without being explicitly coded, while being possible to opt-out.
+3. Provides some utility functions for further building field-specific libraries.
+
+See the documentation online for details on how to use this package as well as examples highlighting its usefulness.
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..ed92fe2
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
+DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..a9cdcfe
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,25 @@
+cd(@__DIR__)
+
+using ProcessBasedModelling
+
+import Downloads
+Downloads.download(
+    "https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/build_docs_with_style.jl",
+    joinpath(@__DIR__, "build_docs_with_style.jl")
+)
+include("build_docs_with_style.jl")
+
+pages =  [
+    "Introduction" => "index.md",
+    "Overarching tutorial" => "tutorial.md",
+    "Contents" => "contents.md",
+    "Animations, GUIs, Visuals" => "visualizations.md",
+    "Contributor Guide" => "contributors_guide.md",
+]
+
+build_docs_with_style(pages, ProcessBasedModelling;
+    authors = "George Datseris <datseris.george@gmail.com>",
+    # We need to remove the cross references because we don't list here
+    # the whole `DynamicalSystem` API...
+    warnonly = [:doctest, :missing_docs, :cross_references],
+)
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..b7fe99b
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,164 @@
+```@docs
+ProcessBasedModelling
+```
+
+!!! note "Basic familiarity with ModelingToolkit.jl"
+    These docs assume that you have some basic familiarity with ModelingToolkit.jl. If you don't going through the introductory tutorial of [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) should be enough to get you started!
+
+
+## Usage
+
+In ProcessBasedModelling.jl, each variable is governed by a "process", which in code is just a symbolic expression from ModellingToolkit.jl.
+To coupled the variable with the process it is governed by, a user either defines simple equations of the form `variable ~ expression`, or creates an instance of [`Process`](@ref) if the left-hand-side of the equation needs to be anything more complex.
+Once all the processes about the physical system are collected, they are given as a `Vector` to the [`processes_to_mtkmodel`](@ref) central function. This function also defines what quantifies as a "process".
+
+## Example
+
+Let's say we want to build the system of equations
+
+```math
+\dot{z} = x^2 - z \\
+\dot{x} = 0.1y \\
+y = z - x
+```
+
+symbolically using ModelingToolkit.jl (**MTK**). We define
+
+```@example MAIN
+using ModelingToolkit
+using OrdinaryDiffEq: Tsit5
+
+@variables z(t) = 0.0
+@variables x(t) # no default value
+@variables y(t) = 0.0
+```
+ProcessBasedModelling.jl (**PBM**) strongly recommends that all defined variables have a default value at definition point. Here we didn't do this for ``x`` to illustrate what how such an "omission" will be treated by **PBM**.
+
+To make the equations we want, we can use MTK directly, and call
+```@example MAIN
+@variables t
+
+eqs = [
+  Differential(t)(z) ~ x^2 - z
+  Differential(x) ~ 0.1y
+  y ~ z - x
+]
+
+model = ODESystem(eqs, t; name = :example)
+
+equations(model)
+```
+
+All good. Now, if we missed the process for one variable (because of our own error/sloppyness/very-large-codebase), MTK will not throw an error at model construction,
+
+```@example MAIN
+model = ODESystem(eqs[1:2], t; name = :example)
+model = structural_simplify(model)
+equations(model)
+```
+
+only at the construction of the "problem" (here the `ODEProblem`)
+
+```@example MAIN
+try
+  prob = ODEProblem(model)
+catch e
+  return e.msg
+end
+```
+
+Interestingly, the error is wrong. ``x`` is defined and has an equation, at least on the basis of our scientific reasoning. However ``y`` that ``x`` introduced does not have an equation. Moreover, in our experience these errors messages become increasingly less useful when a model has many equations and/or variables, as many variables get cited as "missing" from the variable map even when only one should be.
+
+**PBM** resolves these problems and always gives accurate error messages. This is because on top of the variable map that MTK constructs automatically, **PBM** requires the user to implicitly create a map of variables to processes that govern said variables. **PBM** creates the map automatically, the only thing the user has to do is to define the equations in terms of what [`processes_to_mtkmodel`](@ref) wants (which are either [`Process`](@ref)es or `Equation`s as above).
+Here is what the user defines to make the same system of equations:
+
+```@example MAIN
+processes = [
+    ExpRelaxation(z, x^2), # introduces x variable
+    TimeDerivative(x, 0.1*y),   # introduces y variable
+    y ~ z-x,                    # can be an equation because LHS is single variable
+]
+```
+
+Internally, all of these
+
+which is the given to
+```@example MAIN
+model = processes_to_mtkmodel(processes)
+equations(model)
+```
+
+Notice that the resulting **MTK** model is not `structural_simplify`-ed, to allow composing it with other models.
+
+Now, in contrast to before, if we "forgot" a process, **PBM** will react accordingly. For example, we forgot the 2nd process, then the construction will error informatively, telling us exactly which variable is missing, and because of which processes it is missing:
+```@example MAIN
+try
+  model = processes_to_mtkmodel(processes[[1, 3]])
+catch e
+  return e.msg
+end
+```
+
+If instead we "forgot" the ``y`` process, **PBM** will not error, but instead warn, and make ``y`` a named parameter:
+```@example MAIN
+model = processes_to_mtkmodel(processes[1:2])
+equations(model)
+```
+
+```@example MAIN
+parameters(model)
+```
+
+### Special handling of timescales
+
+In dynamical systems modelling the timescale associated with a process is a special parameter. That is why, if a timescale is given for either the [`TimeDerivative`](@ref) or [`ExpRelaxation`](@ref) processes, it is converted to a named `@parameter` by default:
+
+```@example MAIN
+processes = [
+    ExpRelaxation(z, x^2, 2.0),  # third argument is the timescale
+    TimeDerivative(x, 0.1*y, 0.5),
+    y ~ z-x,
+]
+
+model = processes_to_mtkmodel(processes)
+equations(model)
+```
+
+```@example MAIN
+parameters(model)
+```
+
+This special handling is also why each process explicitly declares a timescale via the [`timescale`](@ref) function that one can optionally extend.
+
+
+## Main API function
+
+```@docs
+processes_to_mtkmodel
+```
+
+## [Predefined `Process` subtypes](@id predefined_processes)
+
+```@docs
+ParameterProcess
+TimeDerivative
+ExpRelaxation
+```
+
+## `Process` API
+
+This API describes how you can implement your own `Process` subtype, if the [existing predefined subtypes](@ref predefined_processes) don't fit your bill!
+
+```@docs
+Process
+rhs
+timescale
+NoTimeVariability
+```
+
+## Utility functions
+
+```@docs
+has_variable
+new_derived_named_parameter
+```
diff --git a/src/API.jl b/src/API.jl
new file mode 100644
index 0000000..0d2a491
--- /dev/null
+++ b/src/API.jl
@@ -0,0 +1,60 @@
+"""
+A process subtype `p::Process` extends the functions:
+- `lhs_variable(p)` which returns the variable the process describes (left-hand-side variable)
+- `rhs(p)` which is the right-hand-side expression, i.e., the "actual" process.
+- (optional) `timescale`, which defaults to [`NoTimeVariability`](@ref).
+- (optional) `lhs(p)` which returns the left-hand-side. Let `τ = timescale(p)`.
+  Then `lhs(p)` behavior depends on `τ` as follows:
+  - Just `lhs_variable(p)` if `τ == NoTimeVariability()`.
+  - `Differential(t)(p)` if `τ == nothing`.
+  - `τ_var*Differential(t)(p)` if `τ isa Union{Real, Num}`. If real,
+    a new named parameter `τ_var` is created that has the prefix `:τ_` and then the
+    lhs-variable name and has default value `τ`. Else if `Num`, `τ_var = τ` as given.
+"""
+abstract type Process end
+
+"""
+    NoTimeVariability()
+
+Singleton value that is the default output of the [`timescale`](@ref) function
+for variables that do not vary in time autonomously (i.e., no d/dt derivative).
+"""
+struct NoTimeVariability end
+
+function lhs_variable(p::Process)
+    if !hasfield(typeof(p), :variable)
+        error("`lhs_variable` not defined for process $(nameof(typeof(p))).")
+    else
+        return p.variable
+    end
+end
+
+timescale(::Process) = NoTimeVariability()
+
+function lhs(p::Process)
+    τ = timescale(p)
+    v = lhs_variable(p)
+    if isnothing(τ) # time variability exists but timescale is nonexistent (unity)
+        return Differential(t)(v)
+    elseif τ isa NoTimeVariability || iszero(τ) # no time variability
+        return v
+    else # τ is either Num or Real
+        τvar = new_derived_named_parameter(v, τ, "τ", false)
+        return τvar*Differential(t)(v)
+    end
+end
+function rhs(p::Process)
+    error("Right-hand side (`rhs`) is not defined for process $(nameof(typeof(p))).")
+end
+
+rhs(e::Equation) = e.rhs
+lhs(e::Equation) = lhs_variable(e)
+function lhs_variable(e::Equation)
+    x = e.lhs
+    # we first check whether `x` is a variable
+    if !is_variable(x)
+        throw(ArgumentError("In given equation $(e), the left-hand-side does "*
+        "not represent a variable."))
+    end
+    return x
+end
diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl
new file mode 100644
index 0000000..a5dc1fc
--- /dev/null
+++ b/src/ProcessBasedModelling.jl
@@ -0,0 +1,38 @@
+module ProcessBasedModelling
+
+# Use the README as the module docs
+@doc let
+    path = joinpath(dirname(@__DIR__), "README.md")
+    include_dependency(path)
+    read(path, String)
+end ProcessBasedModelling
+
+using Reexport
+@reexport using ModelingToolkit
+
+@variables t # independent variable (time)
+
+include("API.jl")
+include("utils.jl")
+include("make.jl")
+include("processes_basic.jl")
+
+
+# TODO: In MAKE, make it so that if a variable does not have a process
+# a constant process is created for it if it has a default value.
+# Add a keyword `use_default` which would warn if no process but there is default
+# and otherwise would error.
+# TODO: Make an "addition process" that adds to processes
+# It checks whether they target the same variable
+# TODO: Package should compose with ODESystem
+# so that component-based modelling can be utilized as well.
+
+
+# TODO: Perhaps not don't export `t, rhs`?
+export t
+export lhs_variable, rhs, timescale, NoTimeVariability
+export Process, ParameterProcess, TimeDerivative, ExpRelaxation
+export processes_to_mtkmodel
+export new_derived_named_parameter
+
+end
diff --git a/src/make.jl b/src/make.jl
new file mode 100644
index 0000000..1f186fa
--- /dev/null
+++ b/src/make.jl
@@ -0,0 +1,147 @@
+"""
+    processes_to_mtkmodel(processes::Vector, default::Vector = []; kw...)
+
+Construct a ModelingToolkit.jl model using the provided `processes` and `default` processes.
+
+`processes` is a vector whose elements can be:
+
+- Any instance of a subtype of [`Process`](@ref).
+- An `Equation` which is of the form `variable ~ expression` with `variable` a single
+  variable resulting from an `@variables` call.
+- A vector of the above two, which is then expanded. This allows the convenience of
+  functions representing a physical process that may require many equations to be defined.
+
+`default` is a vector that can contain the first two possibilities only
+as it contains default processes that may be assigned to variables introduced in `processes`
+but they don't themselves have an assigned process.
+
+It is expected that downstream packages that use ProcessBasedModelling.jl to make a
+field-specific library implement a 1-argument version of `processes_to_mtkmodel`,
+or provide a wrapper function for it, and add a default value for `default`.
+"""
+function processes_to_mtkmodel(_processes, _default = []; type = ODESystem, name = nameof(type))
+    processes = expand_multi_processes(_processes)
+    default = default_dict(_default)
+    # Setup: obtain lhs-variables so we can track new variables that are not
+    # in this vector
+    lhs_vars = map(lhs_variable, processes)
+    ensure_unique_vars(lhs_vars)
+    # First pass: generate equations from given processes
+    # and collect variables without equations
+    incomplete = Num[]
+    introduced = Dict{Num, Num}()
+    eqs = Equation[]
+    for proc in processes
+        append_incomplete_variables!(incomplete, introduced, lhs_vars, proc)
+        # add the generated equation in the pool of equations
+        push!(eqs, lhs(proc) ~ rhs(proc))
+    end
+    # Second pass: attempt to add default processes to incomplete variables
+    # throw an error if default does not exist
+    while !isempty(incomplete) # using a `while` allows us check variables added by the defaults
+        added_var = popfirst!(incomplete)
+        if haskey(default, added_var)
+            # Then obtain default process
+            def_proc = default[added_var]
+            # add the equation to the equation pool
+            push!(eqs, lhs(def_proc) ~ rhs(def_proc))
+            # Also ensure that potentially new processes introduced by the default process
+            # are taken care of and have a default value
+            push!(lhs_vars, lhs_variable(def_proc))
+            append_incomplete_variables!(incomplete, introduced, lhs_vars, def_proc)
+        else
+            def_val = default_value(added_var) # utilize default value (if possible)
+            if !isnothing(def_val)
+                @warn("""
+                Variable $(added_var) was introduced in process of variable $(introduced[added_var]).
+                However, a process for $(added_var) was not provided,
+                and there is no default process for it either.
+                Since it has a default value, we make it a parameter by adding a process:
+                `ParameterProcess($(added_var))`.
+                """)
+                parproc = ParameterProcess(added_var)
+                push!(eqs, lhs(parproc) ~ rhs(parproc))
+                push!(lhs_vars, added_var)
+            else
+                throw(ArgumentError("""
+                Variable $(added_var) was introduced in process of variable $(introduced[added_var]).
+                However, a process for $(added_var) was not provided,
+                there is no default process for, and it doesn't have a default value.
+                Please provide a process for variable $(added_var).
+                """))
+            end
+        end
+    end
+    sys = type(eqs, t; name)
+    return sys
+end
+# version without given processes
+function processes_to_mtkmodel(; kwargs...)
+    return processes_to_mtkmodel(collect(values(default_processes())); kwargs...)
+end
+
+function expand_multi_processes(procs::Vector)
+    # Expand vectors of processes or ODESystems
+    !any(p -> p isa Vector, procs) && return procs
+    expanded = deepcopy(procs)
+    idxs = findall(p -> p isa Vector, procs)
+    multiprocs = expanded[idxs]
+    deleteat!(expanded, idxs)
+    for mp in multiprocs
+        append!(expanded, mp)
+    end
+    return expanded
+end
+
+function default_dict(processes)
+    default = Dict{Num, Any}()
+    for proc in processes
+        key = lhs_variable(proc)
+        if haskey(default, key)
+            throw(ArgumentError("More than 1 processes for variable $(key) in default processes."))
+        end
+        default[key] = proc
+    end
+    return default
+end
+
+# Add variables to `incomplete` from expression `r`, provided they are not in `lhs_vars`
+# already. Also record which variables added them
+function append_incomplete_variables!(incomplete, introduced, lhs_vars, process)
+    proc_vars = filter(is_variable, get_variables(rhs(process)))
+    newvars = setdiff(proc_vars, lhs_vars)
+    # Store newly introduced variables without a lhs
+    for nv in newvars
+        # Note we can't use `(nv ∈ incomplete)`, that creates a symbolic expression
+        any(isequal(nv), incomplete) && continue # skip if is already recorded
+        numnv = (nv isa Num) ? nv : Num(nv)
+        Symbol(nv) == :t && continue # Time-dependence is not a new variable!
+        push!(incomplete, numnv)
+        # Also record which variable introduced it
+        introduced[numnv] = lhs_variable(process)
+    end
+    return
+end
+
+function ensure_unique_vars(lhs_vars)
+    nonun = nonunique(lhs_vars)
+    isempty(nonun) || error("The following variables have more than one processes assigned to them: $(nonun)")
+    return
+end
+
+function nonunique(x::AbstractArray{T}) where T
+    uniqueset = Set{T}()
+    duplicatedset = Set{T}()
+    duplicatedvector = Vector{T}()
+    for i in x
+        if(i in uniqueset)
+            if !(i in duplicatedset)
+                push!(duplicatedset, i)
+                push!(duplicatedvector, i)
+            end
+        else
+            push!(uniqueset, i)
+        end
+    end
+    duplicatedvector
+end
\ No newline at end of file
diff --git a/src/processes_basic.jl b/src/processes_basic.jl
new file mode 100644
index 0000000..f5d11c5
--- /dev/null
+++ b/src/processes_basic.jl
@@ -0,0 +1,80 @@
+"""
+    ParameterProcess(variable, value = default_value(variable)) <: Process
+
+The simplest process which equates a given `variable` to a constant value
+that is encapsulated in a parameter. If `value isa Real`, then
+hence, a named parameter with the name of `variable` and `_0` appended is created.
+Else, if `valua isa Num` then it is taken as the paremeter directly.
+
+Example:
+```julia
+@variables T(t) = 0.5
+proc = ParameterProcess(T)
+```
+will create the equation `T ~ T_0`, where `T_0` is a `@parameter` with default value `0.5`.
+"""
+struct ParameterProcess <: Process
+    variable
+    value
+    function ParameterProcess(variable, value)
+        var = new_derived_named_parameter(variable, value, "0")
+        return new(variable, var)
+    end
+end
+ParameterProcess(var) = ParameterProcess(var, default_value(var))
+lhs_variable(c::ParameterProcess) = c.variable
+rhs(cv::ParameterProcess) = cv.value
+
+"""
+    TimeDerivative(variable, expression [, τ])
+
+The second simplest process that equates the time derivative of
+the `variable` to the given `expression`
+while providing some conveniences over manually constructing an `Equation`.
+
+It creates the equation `τ_\$(variable) Differential(t)(variable) ~ expression`
+by constructing a new `@parameter` with default value `τ`
+(if `τ` is already a `@parameter`, it is used as-is).
+If `τ` is not given, then 1 is used at its place and no parameter is created.
+
+Note that if `iszero(τ)`, then the process `variable ~ expression` is created.
+"""
+struct TimeDerivative <: Process
+    variable
+    expression
+    timescale
+end
+TimeDerivative(a, b) = TimeDerivative(a, b, nothing)
+timescale(e::TimeDerivative) = e.timescale
+rhs(e::TimeDerivative) = e.expression
+
+
+"""
+    ExpRelaxation(variable, expression [, τ]) <: Process
+
+A common process for creating an exponential relaxation of `variable` towards the
+given `expression`, with timescale `τ`. It creates the equation:
+```
+τn*Differential(t)(variable) ~ expression - variable
+```
+Where `τn` is a new named `@parameter` with the value of `τ`
+and name `τ_(\$(variable))`. If instead `τ` is `nothing`, then 1 is used in its place.
+If `iszero(τ)`, then the equation `variable ~ expression` is created instead.
+
+The convenience function
+```julia
+ExpRelaxation(process, τ)
+```
+allows converting an existing process (or equation) into an exponential relaxation
+by using the `rhs` as the `expression` in the equation above.
+"""
+struct ExpRelaxation <: Process
+    variable
+    expression
+    timescale
+end
+ExpRelaxation(v, e) = ExpRelaxation(v, e, nothing)
+ExpRelaxation(proc::Union{Process,Equation}, τ) = ExpRelaxation(lhs_variable(proc), rhs(proc), τ)
+
+timescale(e::ExpRelaxation) = e.timescale
+rhs(e::ExpRelaxation) = iszero(e.timescale) ? e.expression : e.expression - e.variable
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000..92612bb
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,92 @@
+export print_system_info, has_variable, default_value
+export new_derived_named_parameter, @named_parameters
+
+"""
+    has_variable(eq, var)
+
+Return `true` if variable `var` exists in the equation(s) `eq`, `false` otherwise.
+Function works irrespectively if `var` is an `@variable` or `@parameter`.
+"""
+function has_variable(eq::Equation, var)
+    vars = get_variables(eq)
+    return any(isequal(var), vars)
+end
+has_variable(eqs, var) = any(eq -> has_variable(eq, var), eqs)
+
+default_value(x) = x
+default_value(x::Num) = default_value(x.val)
+function default_value(x::ModelingToolkit.SymbolicUtils.Symbolic)
+    if haskey(x.metadata, ModelingToolkit.Symbolics.VariableDefaultValue)
+        return x.metadata[ModelingToolkit.Symbolics.VariableDefaultValue]
+    else
+        @warn("No default value assigned to variable/parameter $(x).")
+        return nothing
+    end
+end
+# trick to get default values for state variables:
+# Base.Fix1(getindex, ModelingToolkit.defaults(ssys)).(states(ssys))
+# while `defaults` returns all default assignments.
+
+is_variable(x::Num) = is_variable(x.val)
+function is_variable(x)
+    if x isa ModelingToolkit.SymbolicUtils.Symbolic
+        if isnothing(x.metadata)
+            return false
+        end
+        if haskey(x.metadata, ModelingToolkit.Symbolics.VariableSource)
+            src = x.metadata[ModelingToolkit.Symbolics.VariableSource]
+            return first(src) == :variables
+        end
+    end
+    return false
+end
+
+"""
+    new_derived_named_parameter(variable, value, extra::String, suffix = true)
+
+If `value isa Num`, return `value`. Otherwise, create a new MTK `@parameter`
+whose name is created from `variable` by adding the `extra` string.
+If `suffix = true` the extra is added at the end after a `_`. Otherwise
+it is added at the start, then a `_` and then the variable name.
+"""
+new_derived_named_parameter(v, value::Num, extra, suffix = true) = value
+function new_derived_named_parameter(v, value::Real, extra, suffix = true)
+    n = string(ModelingToolkit.getname(v))
+    newstring = if suffix
+        n*"_"*extra
+    else
+        extra*"_"*n
+    end
+    new_derived_named_parameter(newstring, value)
+end
+function new_derived_named_parameter(newstring::String, value::Real)
+    varsymbol = Symbol(newstring)
+    dummy = (@parameters $(varsymbol) = value)
+    return first(dummy)
+end
+
+"""
+    @named_parameters vars...
+
+Convert all Julia variables `vars` into `@parameters` with name the same as `vars`
+and default value the same as the value of `vars`.
+"""
+macro named_parameters(vars...)
+    return quote
+        out = []
+        for v in vars
+            res = (ModelingToolkit.toparam)((Symbolics.wrap)((Symbolics.SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}($(QuoteNode(v))), $(esc(v))), Symbolics.VariableSource, (:parameters, $(QuoteNode(v))))))
+            push!(out, res)
+        end
+        $(out...)
+    end
+end
+
+# This may help:
+
+# julia> @macroexpand @parameters A=0.1 B=0.1
+# quote
+#     A = (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}(:A), 0.1), Symbolics.VariableSource, (:parameters, :A))))
+#     B = (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}(:B), 0.1), Symbolics.VariableSource, (:parameters, :B))))
+#     [A, B]
+# end
\ No newline at end of file
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..d170eec
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,6 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
new file mode 100644
index 0000000..a2d025b
--- /dev/null
+++ b/test/runtests.jl
@@ -0,0 +1,95 @@
+using ProcessBasedModelling
+using Test
+using OrdinaryDiffEq
+
+@testset "construction + evolution" begin
+    # The model, as defined below, is bistable due to ice albedo feedback
+    # so two initial conditions should go to two attractors
+    # If that's the case, we are sure model construction was valid
+
+    # First, make some default processes
+    @variables T(t) = 300.0       # temperature, in Kelvin
+    @variables α(t) = 0.3         # albedo of ice, unitless
+    @variables ε(t) = 0.5         # effective emissivity, unitless
+    solar_constant = 340.25 # W/m^2, already divided by 4
+    σ_Stefan_Boltzman = 5.670374419e-8    # stefan boltzman constant
+
+    Base.@kwdef struct HeatBalance <: Process
+        c_T = 5e8
+    end
+    ProcessBasedModelling.lhs_variable(::HeatBalance) = T
+    ProcessBasedModelling.timescale(proc::HeatBalance) = proc.c_T/solar_constant
+    function ProcessBasedModelling.rhs(::HeatBalance)
+        absorbed_shortwave = 1 - α
+        emitted_longwave = ε*(σ_Stefan_Boltzman/solar_constant)*T^4
+        return absorbed_shortwave - emitted_longwave
+    end
+
+    processes = [
+        TanhProcess(α, T, 0.7, 0.289, 10.0, 274.5),
+        TanhProcess(ε, T, 0.5, 0.41, 2.0, 288.0),
+        HeatBalance()
+    ]
+
+    sys = processes_to_mtkmodel(processes)
+    @test sys isa ODESystem
+    @test length(states(sys)) == 3
+
+    sys = structural_simplify(sys)
+    @test length(states(sys)) == 1
+
+    u0s = [[300.0], [100.0]]
+    ufs = []
+    for u0 in u0s
+        p = ODEProblem(sys, u0, (0.0, 1000.0*365*24*60*60.0))
+        sol = solve(p, Tsit5())
+        push!(ufs, sol.u[end])
+    end
+
+    @test ufs[1] ≈ [319] atol = 1
+    @test ufs[2] ≈ [245] atol = 1
+end
+
+@testset "add missing processes" begin
+    @variables z(t) = 0.0
+    @variables x(t) # no default value
+    @variables y(t) = 0.0
+
+    procs = [
+        ExpRelaxation(z, x^2, 1.0), # introduces x and y variables
+        TimeDerivative(x, 0.1*y),   # introduces y variable!
+        y ~ z-x,                    # is an equation, not a process!
+    ]
+
+    @testset "only one variable and no defaults" begin
+        @test_throws ArgumentError processes_to_mtkmodel(procs[1:1])
+    end
+
+    @testset "first two, still missing y, but it has default" begin
+        model = @test_logs (:warn, r"\W*((?i)Variable(?-i))\W*") processes_to_mtkmodel(procs[1:2])
+        @test length(states(model)) == 3
+    end
+
+    @testset "first with default the third; missing x" begin
+        @test_throws ArgumentError processes_to_mtkmodel(procs[1:1], procs[3:3])
+    end
+
+    @testset "first with default the second; y gets contant value a different warning" begin
+        model = @test_logs (:warn, r"\W*((?i)parameter(?-i))\W*") processes_to_mtkmodel(procs[1:1], procs[2:2])
+        @test length(states(model)) == 3
+    end
+
+    @testset "all three processes given" begin
+        sys = processes_to_mtkmodel(procs[1:1], procs[2:3])
+        @test length(states(sys)) == 3
+        sys = processes_to_mtkmodel(procs[1:2], procs[3:3])
+        @test length(states(sys)) == 3
+        sys = processes_to_mtkmodel(procs[1:3])
+        @test length(states(sys)) == 3
+        @test length(states(structural_simplify(sys))) == 2
+    end
+end
+
+@testset "extending default processes" begin
+    # API not yet finished on this one
+end