diff --git a/README.md b/README.md
index 252e8a7..f928542 100644
--- a/README.md
+++ b/README.md
@@ -52,7 +52,7 @@ Here `BNGNetwork` is a type specifying the file format that is being loaded.
 `prnbng` is a `ParsedReactionNetwork` structure with the following fields:
 
   - `rn`, a Catalyst `ReactionSystem`
-  - `u₀`, a `Dict` mapping initial condition symbolic variables to numeric values
+  - `u0`, a `Dict` mapping initial condition symbolic variables to numeric values
     and/or symbolic expressions.
   - `p`, a `Dict` mapping parameter symbolic variables to numeric values and/or
     symbolic expressions.
@@ -221,4 +221,4 @@ initialconditionf = "PATH/TO/FILE"
 networkf = "PATH/TO/FILE"
 rssarn = loadrxnetwork(RSSANetwork(), "RSSARxSys", initialconditionf, networkf)
 ```
-Here `RSSANetwork` specifies the type of the file to parse, and `RSSARxSys` gives the type of the generated `reaction_network`. `rssarn` is again a `ParsedReactionNetwork`, but only the `rn` and `u₀` fields will now be relevant (the remaining fields will be set to `nothing`). -->
+Here `RSSANetwork` specifies the type of the file to parse, and `RSSARxSys` gives the type of the generated `reaction_network`. `rssarn` is again a `ParsedReactionNetwork`, but only the `rn` and `u0` fields will now be relevant (the remaining fields will be set to `nothing`). -->
diff --git a/docs/src/index.md b/docs/src/index.md
index 2452885..af3106c 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -64,7 +64,7 @@ Here `BNGNetwork` is a type specifying the file format that is being loaded.
 `prnbng` is a `ParsedReactionNetwork` structure with the following fields:
 
   - `rn`, a Catalyst `ReactionSystem`
-  - `u₀`, a `Dict` mapping initial condition symbolic variables to numeric values
+  - `u0`, a `Dict` mapping initial condition symbolic variables to numeric values
     and/or symbolic expressions.
   - `p`, a `Dict` mapping parameter symbolic variables to numeric values and/or
     symbolic expressions.
diff --git a/examples/not_updated_yet/bcr_odes_sparse.jl b/examples/not_updated_yet/bcr_odes_sparse.jl
index 6caf43a..950c681 100644
--- a/examples/not_updated_yet/bcr_odes_sparse.jl
+++ b/examples/not_updated_yet/bcr_odes_sparse.jl
@@ -32,7 +32,7 @@ reset_timer!(to)
 # BioNetGen network
 @timeit to "bionetgen" prnbng=loadrxnetwork(BNGNetwork(), string(networkname, "bng"), fname);
 rnbng = prnbng.rn;
-u0 = prnbng.u₀;
+u0 = prnbng.u0;
 p = prnbng.p;
 @timeit to "baddodes" addodes!(rnbng; build_jac = build_jac, zeroout_jac = zeroout_jac,
                                sparse_jac = sparse_jac, build_symfuncs = false,
diff --git a/examples/not_updated_yet/test_network.jl b/examples/not_updated_yet/test_network.jl
index 92f2a33..100ba79 100644
--- a/examples/not_updated_yet/test_network.jl
+++ b/examples/not_updated_yet/test_network.jl
@@ -21,18 +21,18 @@ reset_timer!(to)
 @timeit to "netgen" prn=loadrxnetwork(RSSANetwork(), networkname, speciesf, rxsf;
                                       printrxs = false)
 rn = prn.rn
-initialpop = prn.u₀
+initialpop = prn.u0
 println("network parsed")
 
 # get the BioNetGen reaction network
 @timeit to "bionetgen" prnbng=loadrxnetwork(BNGNetwork(), string(networkname, "bng"),
                                             bngfname)
 rnbng = prnbng.rn;
-u₀ = prnbng.u₀;
+u0 = prnbng.u0;
 p = prnbng.p;
 shortsymstosyms = prnbng.symstonames;
 
-u0 = round.(Int, u₀)
+u0 = round.(Int, u0)
 println(typeof(u0))
 
 # one simulation
diff --git a/examples/not_updated_yet/test_odes.jl b/examples/not_updated_yet/test_odes.jl
index a7e30bd..ad25ca2 100644
--- a/examples/not_updated_yet/test_odes.jl
+++ b/examples/not_updated_yet/test_odes.jl
@@ -16,7 +16,7 @@ reset_timer!(to)
 # get the reaction network
 @timeit to "netgen" prnbng=loadrxnetwork(BNGNetwork(), string(networkname, "bng"), fname);
 rn = prnbng.rn;
-u0 = prnbng.u₀;
+u0 = prnbng.u0;
 p = prnbng.p;
 @timeit to "addodes" addodes!(rn; build_jac = build_jac, sparse_jac = sparse_jac,
                               build_symfuncs = false, build_paramjac = false)
diff --git a/examples/not_updated_yet/test_odes_rssa.jl b/examples/not_updated_yet/test_odes_rssa.jl
index 47397b2..601c295 100644
--- a/examples/not_updated_yet/test_odes_rssa.jl
+++ b/examples/not_updated_yet/test_odes_rssa.jl
@@ -20,7 +20,7 @@ reset_timer!(to)
 @timeit to "netgen" prn=loadrxnetwork(RSSAFile(), networkname, speciesf, rxsf;
                                       printrxs = false)
 rn = prn.rn;
-initialpop = prn.u₀;
+initialpop = prn.u0;
 @timeit to "addodes" addodes!(rn; build_jac = false, build_symfuncs = false,
                               build_paramjac = false)
 @timeit to "ODEProb" oprob=ODEProblem(rn, convert.(Float64, initialpop), (0.0, tf))
diff --git a/src/ReactionNetworkImporters.jl b/src/ReactionNetworkImporters.jl
index 645e97f..4de6c94 100644
--- a/src/ReactionNetworkImporters.jl
+++ b/src/ReactionNetworkImporters.jl
@@ -24,7 +24,7 @@ struct ParsedReactionNetwork
     rn::ReactionSystem
 
     "Dict mapping initial condition symbolic variables to values."
-    u₀::Any
+    u0::Any
 
     "Dict mapping parameter symbolic variables to values."
     p::Any
@@ -35,9 +35,9 @@ struct ParsedReactionNetwork
     "Dict from group name (as string) to corresponding symbolic variable"
     groupstosyms::Any
 end
-function ParsedReactionNetwork(rn::ReactionSystem; u₀ = nothing, p = nothing,
+function ParsedReactionNetwork(rn::ReactionSystem; u0 = nothing, p = nothing,
                                varstonames = nothing, groupstosyms = nothing)
-    ParsedReactionNetwork(rn, u₀, p, varstonames, groupstosyms)
+    ParsedReactionNetwork(rn, u0, p, varstonames, groupstosyms)
 end
 
 export BNGNetwork, MatrixNetwork, ParsedReactionNetwork, ComplexMatrixNetwork
@@ -48,4 +48,27 @@ include("parsing_routines_matrixnetworks.jl")
 
 export loadrxnetwork
 
+
+# Overload ensuring that u0 and u₀ can be used interchangeably.
+# (introduced when the u₀ field was changed to u0)
+# Should be deleted whenever u₀ is fully deprecated.
+
+# Ensures that `prnbng.u₀` works.
+function Base.getproperty(prnbng::ParsedReactionNetwork, name::Symbol)
+    if name === :u₀
+        return getfield(prnbng, :u0)
+    else
+        return getfield(prnbng, name)
+    end
+end
+
+# Ensures that `prnbng.u₀ = ...` works.
+function Base.setproperty!(prnbng::ParsedReactionNetwork, name::Symbol, x)
+    if name === :u₀
+        return setfield!(prnbng, :u0, x)
+    else
+        return setfield!(prnbng, name, x)
+    end
+end
+
 end # module
diff --git a/src/parsing_routines_bngnetworkfiles.jl b/src/parsing_routines_bngnetworkfiles.jl
index d31e778..dce131e 100644
--- a/src/parsing_routines_bngnetworkfiles.jl
+++ b/src/parsing_routines_bngnetworkfiles.jl
@@ -279,10 +279,10 @@ function loadrxnetwork(ft::BNGNetwork, rxfilename; name = gensym(:ReactionSystem
     rn = ReactionSystem(rxs, t, specs, ps; name = name, observed = obseqs,
                         defaults = defmap, kwargs...)
 
-    # get numeric values for parameters and u₀
+    # get numeric values for parameters and u0
     sm = speciesmap(rn)
     @assert all(sm[funcsym(sym, t)] == i for (i, sym) in enumerate(idstoshortsyms))
 
-    ParsedReactionNetwork(rn; u₀ = u0map, p = pmap, varstonames = shortsymstosyms,
+    ParsedReactionNetwork(rn; u0 = u0map, p = pmap, varstonames = shortsymstosyms,
                           groupstosyms = groupstosyms)
 end
diff --git a/test/test_higherorder_odes.jl b/test/test_higherorder_odes.jl
index f86d98d..3c7c832 100644
--- a/test/test_higherorder_odes.jl
+++ b/test/test_higherorder_odes.jl
@@ -21,6 +21,9 @@ prnbng = loadrxnetwork(BNGNetwork(), fname)
 rn = prnbng.rn
 boprob = ODEProblem(rn, Float64[], (0.0, tf), Float64[])
 
+# Test that u0 == u₀ (used when the u0 indexing was introduced).
+@test isequal(prnbng.u0, prnbng.u₀)
+
 # BNG simulation data testing
 @unpack A = rn
 Asol = gdatdf[!, :A]
diff --git a/test/test_nullrxs_odes.jl b/test/test_nullrxs_odes.jl
index 332251b..6af8c54 100644
--- a/test/test_nullrxs_odes.jl
+++ b/test/test_nullrxs_odes.jl
@@ -24,6 +24,9 @@ u0 = Float64[];
 p = Float64[];
 boprob = ODEProblem(rn, u0, (0.0, tf), p)
 
+# Test that u0 == u₀ (used when the u0 indexing was introduced).
+@test isequal(prnbng.u0, prnbng.u₀)
+
 # note solvers run _much_ faster the second time 
 bsol = solve(boprob, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = tf / nsteps);
 
diff --git a/test/test_repressilator_odes.jl b/test/test_repressilator_odes.jl
index bf7a4da..be9e4b2 100644
--- a/test/test_repressilator_odes.jl
+++ b/test/test_repressilator_odes.jl
@@ -23,6 +23,9 @@ u0 = Float64[];
 p = Float64[];
 boprob = ODEProblem(rn, u0, (0.0, tf), p)
 
+# Test that u0 == u₀ (used when the u0 indexing was introduced).
+@test isequal(prnbng.u0, prnbng.u₀)
+
 # BNG simulation data testing
 bngsol = gdatdf[!, :pTetR]