diff --git a/reacnetgenerator/_draw.py b/reacnetgenerator/_draw.py
index 5054e2e94..e13e50e07 100644
--- a/reacnetgenerator/_draw.py
+++ b/reacnetgenerator/_draw.py
@@ -18,6 +18,7 @@
 
 import logging
 import math
+import traceback
 from io import StringIO
 
 import matplotlib.pyplot as plt
@@ -45,10 +46,18 @@ def __init__(self, rng):
         self.pos = rng.pos
         self.nolabel = rng.nolabel
         self.showid = rng.showid
+        self._split = rng.split
 
     def draw(self):
         """Draw the network."""
-        table, name = self._readtable()
+        self._draw()
+        if self._split > 1:
+            for st in range(self._split):
+                self._draw(timeaxis=st)
+
+    def _draw(self, timeaxis=None):
+        table, name = self._readtable(
+            self.tablefilename if timeaxis is None else f"{self.tablefilename}.{timeaxis}")
         species, showname = self._handlespecies(name)
 
         G = nx.DiGraph()
@@ -67,35 +76,35 @@ def draw(self):
                                   tableij)])
         weights = np.array([math.log(G[u][v]['weight']+1)
                             for u, v in G.edges()])
-        widths = [
-            weight / max(weights) * self.widthcoefficient * 2
-            if weight > max(weights) * 0.7 else weight / max(weights) * self.widthcoefficient *
-            0.5 for weight in weights]
-        colors = [
-            self.start_color + weight / max(weights) *
-            (self.end_color - self.start_color) for weight in weights]
+        widths = weights / np.max(weights) * self.widthcoefficient * np.array(
+            [0.5, 2])[(weights > np.max(weights) * 0.7)+0] if weights.size else np.zeros(0)
+        colors = self.start_color + weights[:, np.newaxis] / np.max(weights) * (
+            self.end_color - self.start_color) if weights.size else np.zeros(0)
         try:
-            self.pos = (nx.spring_layout(G) if not self.pos else nx.spring_layout(G, pos=self.pos, fixed=[p for p in self.pos])) if not self.k else (
-                nx.spring_layout(G, k=self.k) if not self.pos else nx.spring_layout(G, pos=self.pos, fixed=[p for p in self.pos], k=self.k))
-            if self.pos:
+            pos = nx.spring_layout(G,
+                                   pos=self.pos if self.pos else None,
+                                   fixed=list(self.pos) if self.pos else None,
+                                   k=self.k)
+            if pos:
                 logging.info("The position of the species in the network is:")
-                logging.info(self.pos)
+                logging.info(pos)
             for with_labels in ([True] if not self.nolabel else [True, False]):
                 nx.draw(
-                    G, pos=self.pos, width=widths, node_size=self.node_size,
+                    G, pos=pos, width=widths, node_size=self.node_size,
                     font_size=self.font_size, with_labels=with_labels,
-                    edge_color=colors, node_color=[self.node_color]*len(self.pos))
+                    edge_color=colors, node_color=[self.node_color]*len(pos))
                 imagefilename = "".join(
                     (("" if with_labels else "nolabel_"), self.imagefilename))
-                with StringIO() as stringio, open(imagefilename, 'w') as f:
+                with StringIO() as stringio, open(imagefilename if timeaxis is None else f"{imagefilename}.{timeaxis}", 'w') as f:
                     plt.savefig(stringio, format='svg')
                     f.write(scour.scour.scourString(stringio.getvalue()))
                 plt.close()
         except Exception as e:
             logging.error(f"Error: cannot draw images. Details: {e}")
+            traceback.print_tb(e.__traceback__)
 
-    def _readtable(self):
-        df = pd.read_csv(self.tablefilename, sep=' ', index_col=0, header=0)
+    def _readtable(self, filename):
+        df = pd.read_csv(filename, sep=' ', index_col=0, header=0)
         return df.values, df.index
 
     def _handlespecies(self, name):
diff --git a/reacnetgenerator/_matrix.py b/reacnetgenerator/_matrix.py
index 26bb7ab46..5ef788730 100644
--- a/reacnetgenerator/_matrix.py
+++ b/reacnetgenerator/_matrix.py
@@ -35,13 +35,14 @@ def generate(self):
             R=[a_ij ], i=1,2,…,100;j=1,2,…,100
         where aij is the number of reactions from species si to sj.
         """
-        allroute = self._getallroute(self.allmoleculeroute)
-        self._printtable(allroute)
+        self._printtable(self._getallroute(self.allmoleculeroute))
+        if self.rng.splitmoleculeroute is not None:
+            for i, smr in enumerate(self.rng.splitmoleculeroute):
+                self._printtable(self._getallroute(smr), timeaxis=i)
         if self.needprintspecies:
             self._printspecies()
 
     def _getallroute(self, allmoleculeroute):
-        allroute = Counter()
         names = self._mname[allmoleculeroute-1]
         names = names[names[:, 0] != names[:, 1]]
         if names.size > 0:
@@ -49,7 +50,7 @@ def _getallroute(self, allmoleculeroute):
             return zip(equations[0].tolist(), equations[1].tolist())
         return []
 
-    def _printtable(self, allroute, maxsize=100):
+    def _printtable(self, allroute, maxsize=100, timeaxis=None):
         species = []
         sortedreactions = sorted(
             allroute, key=lambda d: d[1], reverse=True)
@@ -76,7 +77,7 @@ def _printtable(self, allroute, maxsize=100):
 
         table = np.zeros((maxsize, maxsize), dtype=np.int)
         reactionnumber = np.zeros((2), dtype=np.int)
-        with open(self.reactionfilename, 'w') as f:
+        with open(self.reactionfilename if timeaxis is None else f"{self.reactionfilename}.{timeaxis}", 'w') as f:
             for reaction, n_reaction in sortedreactions:
                 f.write(f"{n_reaction} {'->'.join(reaction)}\n")
                 for i, spec in enumerate(reaction):
@@ -93,7 +94,7 @@ def _printtable(self, allroute, maxsize=100):
 
         df = pd.DataFrame(table[:len(species), :len(
             species)], index=species, columns=species)
-        df.to_csv(self.tablefilename, sep=' ')
+        df.to_csv(self.tablefilename if timeaxis is None else f"{self.tablefilename}.{timeaxis}", sep=' ')
 
     def _searchspecies(self, originspec, sortedreactions, species):
         searchedspecies = []
diff --git a/reacnetgenerator/_path.py b/reacnetgenerator/_path.py
index efbb735ca..fc42220f2 100644
--- a/reacnetgenerator/_path.py
+++ b/reacnetgenerator/_path.py
@@ -46,6 +46,7 @@ def __init__(self, rng):
         self._decompress = rng.decompress
         self._bytestolist = rng.bytestolist
         self._produce = rng.produce
+        self._split = rng.split
         self._mname = None
         self.atomnames = None
 
@@ -67,6 +68,9 @@ def collect(self):
         self._printmoleculename()
         atomeach = self._getatomeach()
         self.rng.allmoleculeroute = self._printatomroute(atomeach)
+        if self._split > 1:
+            splittime = np.array_split(np.arange(self._step), self._split)
+            self.rng.splitmoleculeroute = list([self._printatomroute(atomeach[:, st], timeaxis=i) for i, st in enumerate(splittime)])
         self.rng.mname = self._mname
 
     @abstractmethod
@@ -96,14 +100,14 @@ def _getatomroute(self, item):
             [f"Atom {i} {self.atomname[atomtypei]}: ", " -> ".join(names)])
         return moleculeroute, routestr
 
-    def _printatomroute(self, atomeach):
-        with open(self.atomroutefilename, 'w') as f, Pool(self.nproc, maxtasksperchild=1000) as pool:
+    def _printatomroute(self, atomeach, timeaxis=None):
+        with open(self.atomroutefilename if timeaxis is None else f"{self.atomroutefilename}.{timeaxis}", 'w') as f, Pool(self.nproc, maxtasksperchild=1000) as pool:
             allmoleculeroute = []
             semaphore = Semaphore(self.nproc*150)
             results = pool.imap(self._getatomroute, self._produce(
                 semaphore, enumerate(zip(atomeach, self.atomtype), start=1), ()), 100)
             for moleculeroute, routestr in tqdm(
-                    results, total=self._N, desc="Collect reaction paths",
+                    results, total=self._N, desc="Collect reaction paths" if timeaxis is None else f"Collect reaction paths {timeaxis}",
                     unit="atom"):
                 f.write("".join([routestr, '\n']))
                 if moleculeroute.size > 0:
diff --git a/reacnetgenerator/_reachtml.py b/reacnetgenerator/_reachtml.py
index d9c8d0703..633577c69 100644
--- a/reacnetgenerator/_reachtml.py
+++ b/reacnetgenerator/_reachtml.py
@@ -27,6 +27,7 @@ def __init__(self, rng):
         self._resultfile = rng.resultfilename
         self._imagefile = rng.imagefilename
         self._nproc = rng.nproc
+        self._split = rng.split
         self._templatedict = {
             "speciesshownum": 30,
             "reactionsshownum": 20,
@@ -51,15 +52,15 @@ def _re(cls, smi):
             "C", "[C]").replace(
             "[HH]", "[H]")
 
-    def _readreaction(self):
+    def _readreaction(self, timeaxis=None):
         reaction = []
-        with open(self._reactionfile) as f:
+        with open(self._reactionfile if timeaxis is None else f"{self._reactionfile}.{timeaxis}") as f:
             for line in f:
                 sx = line.split()
                 s = sx[1].split("->")
                 left, right, num = self._re(s[0]), self._re(s[1]), int(sx[0])
                 reaction.append((left, right, num))
-                if len(self._linkreac[left]) < 5:
+                if timeaxis is None and len(self._linkreac[left]) < 5:
                     self._linkreac[left].append(right)
         return reaction
 
@@ -80,29 +81,38 @@ def _convertsvg(cls, smiles):
         svgdata = re.sub(r"""<title>.*?<\/title>""", '', svgdata)
         return smiles, svgdata
 
-    def _readspecies(self):
+    def _readspecies(self, reaction, timeaxis=None):
         specs = []
-        for reac in self._reaction:
+        for reac in reaction:
             for spec in reac[:2]:
                 if spec not in specs:
                     specs.append(spec)
-        with Pool(self._nproc) as pool:
-            results = pool.imap_unordered(self._convertsvg, specs)
-            for spec, svgfile in results:
-                self._svgfiles[spec] = svgfile
-        pool.join()
-        pool.close()
+        if timeaxis is None:
+            with Pool(self._nproc) as pool:
+                results = pool.imap_unordered(self._convertsvg, specs)
+                for spec, svgfile in results:
+                    self._svgfiles[spec] = svgfile
+            pool.join()
+            pool.close()
         return specs
 
     def _readdata(self):
-        self._reaction = self._readreaction()
-        self._specs = self._readspecies()
+        self._reaction = [self._readreaction()]
+        self._specs = [self._readspecies(self._reaction[0])]
+        if self._split > 1:
+            for i in range(self._split):
+                reaction = self._readreaction(timeaxis=i)
+                self._reaction.append(reaction)
+                self._specs.append(self._readspecies(reaction, timeaxis=i))
 
     def _generateresult(self):
-        self._generatenetwork()
+        self._templatedict["network_time"] = [self._generatenetwork()]
+        if self._split > 1:
+            for i in range(self._split):
+                self._templatedict["network_time"].append(self._generatenetwork(timeaxis=i))
         self._generatesvg()
-        self._templatedict["species"] = self._specs
-        self._templatedict["reactions"] = self._reaction
+        self._templatedict["speciestime"] = self._specs
+        self._templatedict["reactionstime"] = self._reaction
         self._templatedict["javascript"] = pkg_resources.resource_string(
             __name__, 'static/webpack/bundle.js').decode()
         self._templatedict["linkreac"] = json.dumps(
@@ -113,16 +123,16 @@ def _generateresult(self):
         with open(self._resultfile, 'w', encoding="utf-8") as f:
             f.write(htmlmin.minify(webpage))
 
-    def _generatenetwork(self):
-        with open(self._imagefile) as f:
+    def _generatenetwork(self, timeaxis=None):
+        with open(self._imagefile if timeaxis is None else f"{self._imagefile}.{timeaxis}") as f:
             svgdata = f.read().strip()
             svgdata = re.sub(r"\d+(\.\d+)?pt", "100%", svgdata, count=2)
             svgdata = re.sub(
                 r"""<(\?xml|\!DOCTYPE|\!\-\-)("[^"]*"|'[^']*'|[^'">])*>""", '',
                 svgdata)
-        self._templatedict["network_svg"] = svgdata
+        return svgdata
 
     def _generatesvg(self):
         self._templatedict["speciessvg"] = list(
             [{"name": spec, "svg": self._svgfiles[spec]}
-             for spec in self._specs])
+             for spec in self._specs[0]])
diff --git a/reacnetgenerator/reacnetgen.py b/reacnetgenerator/reacnetgen.py
index 4527860a0..e2a95cf48 100644
--- a/reacnetgenerator/reacnetgen.py
+++ b/reacnetgenerator/reacnetgen.py
@@ -78,7 +78,7 @@ def __init__(
             needprintspecies=True, speciesfilter=None, node_color=None,
             pos=None, printfiltersignal=False, showid=True, k=None,
             start_color=None, end_color=None, nproc=None, speciescenter=None,
-            n_searchspecies=2):
+            n_searchspecies=2, split=1):
         """Init ReacNetGenerator."""
         print(__doc__)
         print(
@@ -120,6 +120,7 @@ def __init__(
         self.nproc = self._setparam(nproc, cpu_count())
         self.speciescenter = speciescenter
         self.n_searchspecies = n_searchspecies
+        self.split = split
         # define attribute
         self.atomtype = None
         self.step = None
@@ -133,6 +134,7 @@ def __init__(
         self.moleculetempfilename = None
         self.moleculetemp2filename = None
         self.allmoleculeroute = None
+        self.splitmoleculeroute = None
 
     def runanddraw(self, run=True, draw=True, report=True):
         """Analyze the trajectory from MD simulation."""
@@ -280,11 +282,13 @@ def _commandline():
         help='Select atoms in the reaction network, e.g. C', nargs='*')
     parser.add_argument(
         '--stepinterval', help='Step interval', type=int, default=1)
+    parser.add_argument('--split', help='Split number for the time axis', type=int, default=1)
     args = parser.parse_args()
     ReacNetGenerator(
         inputfilename=args.inputfilename, atomname=args.atomname,
         runHMM=not args.nohmm,
         inputfiletype=('lammpsdumpfile' if args.dump else 'lammpsbondfile'),
         nproc=args.nproc, selectatoms=args.selectatoms,
-        stepinterval=args.stepinterval
+        stepinterval=args.stepinterval,
+        split=args.split,
         ).runanddraw()
diff --git a/reacnetgenerator/static/template.html b/reacnetgenerator/static/template.html
index 50adc8b32..3ce2a0d50 100644
--- a/reacnetgenerator/static/template.html
+++ b/reacnetgenerator/static/template.html
@@ -15,8 +15,8 @@
 <body id="page-top">
     <svg class="d-none" id="svgdefs">
         <defs>
-            <path id="narrow" style="text-align:start;line-height:100%;-inkscape-font-specification:Arial Unicode MS" d="M24.35 7.613c-3.035 1.11-5.407 2.908-7.113 5.395h-1.299c.585-1.743 1.567-3.39 2.945-4.938H.649V6.278h18.234c-1.378-1.548-2.36-3.2-2.945-4.956h1.299c1.706 2.487 4.078 4.285 7.114 5.395v.896z"
-                font-size="39.506" font-weight="400" font-family="Arial Unicode MS" />
+            <path id="narrow" style="text-align:start;line-height:100%;" d="M24.35 7.613c-3.035 1.11-5.407 2.908-7.113 5.395h-1.299c.585-1.743 1.567-3.39 2.945-4.938H.649V6.278h18.234c-1.378-1.548-2.36-3.2-2.945-4.956h1.299c1.706 2.487 4.078 4.285 7.114 5.395v.896z"
+                font-size="39.506" font-weight="400"/>
             <circle id="circle" cx="50" cy="50" r="45" stroke="#00f" stroke-width="2" fill="#fff" />
         </defs>
     </svg>
@@ -30,6 +30,21 @@
             </button>
             <div class="collapse navbar-collapse" id="navbarResponsive">
                 <ul class="navbar-nav ml-auto">
+                    {% if network_time|length > 1 %}
+                    <li class="nav-item">
+                        <select id="timeselect">
+                            {% for _ in network_time %}
+                            <option value="{{ loop.index }}">
+                                {% if loop.index > 1 %}
+                                    Time {{ loop.index-1 }}
+                                {% else %}
+                                    All
+                                {% endif %}
+                            </option>
+                            {% endfor %}
+                        </select>
+                    </li>
+                    {% endif %}
                     {% for sec in sections %}
                     <li class="nav-item"><a class="nav-link js-scroll-trigger" href="#{{ sec|lower }}">{{ sec }}</a></li>
                     {% endfor %}
@@ -60,7 +75,8 @@ <h1 class="text-uppercase">Analysis Report</h1>
         </div>
     </header>
 
-    <section id="network" class='page-section bg-white mx-auto text-center'>
+    {% for network_svg in network_time %}
+    <section id="network" class='page-section bg-white mx-auto text-center time time-{{ loop.index }}'>
         <div class="container my-auto">
             <div class="row">
                 <div class="mx-auto w-100 text-center">
@@ -70,6 +86,7 @@ <h2>Reaction Network</h2>
             </div>
         </div>
     </section>
+    {% endfor %}
 
     <svg class="d-none" id="svgspecs">
         <defs>
@@ -88,7 +105,8 @@ <h2>Reaction Network</h2>
         <a href="javascript:savesvg();" class="savesvg">Save</a>
     </div>
 
-    <section id="species" class='page-section bg-white mx-auto'>
+    {% for species in speciestime %}
+    <section id="species" class='page-section bg-white mx-auto time time-{{ loop.index }}'>
         <div class="container">
             <div class="row">
                 <div class="mx-auto text-center">
@@ -114,8 +132,10 @@ <h2 class='text-center'>Species</h2>
                 {% endif %}
             </div>
     </section>
+    {% endfor %}
 
-    <section id="reactions" class='page-section bg-white mx-auto'>
+    {% for reactions in reactionstime %}
+    <section id="reactions" class='page-section bg-white mx-auto time time-{{ loop.index }}'>
         <div class="container">
             <div class="row">
                 <div class="mx-auto text-center">
@@ -149,6 +169,7 @@ <h4 class='text-center'>(sorted by frequency)</h4>
                     {% endif %}
                 </div>
     </section>
+    {% endfor %}
 
     <section class="page-section bg-dark text-white" id="foot">
         <div class="container">
diff --git a/reacnetgenerator/static/webpack/reacnetgen.js b/reacnetgenerator/static/webpack/reacnetgen.js
index 8f2699168..d12f5c571 100644
--- a/reacnetgenerator/static/webpack/reacnetgen.js
+++ b/reacnetgenerator/static/webpack/reacnetgen.js
@@ -62,6 +62,12 @@ $(function () {
         "type": "inline",
         "preloader": false,
     });
+    $("select#timeselect").change(function(){
+        $(".time").hide();
+        $(".time-"+$(this).val()).show();
+    });
+    $(".time").hide();
+    $(".time-1").show();
 });
 
 /**
diff --git a/reacnetgenerator/test/test.json b/reacnetgenerator/test/test.json
index b8d5383bb..f7ac81569 100644
--- a/reacnetgenerator/test/test.json
+++ b/reacnetgenerator/test/test.json
@@ -12,7 +12,8 @@
         ],
         "inputfiletype": "lammpsbondfile",
         "hmm": true,
-        "smiles": true
+        "smiles": true,
+        "split": 10
     },
     {
         "url": [
diff --git a/reacnetgenerator/test/test_reacnetgen.py b/reacnetgenerator/test/test_reacnetgen.py
index 80d4bcdf5..1ae7dd9c6 100644
--- a/reacnetgenerator/test/test_reacnetgen.py
+++ b/reacnetgenerator/test/test_reacnetgen.py
@@ -38,7 +38,9 @@ def reacnetgen(self, request):
             inputfiletype=testparm['inputfiletype'],
             runHMM=testparm['hmm'],
             speciescenter=testparm['speciescenter']
-            if 'speciescenter' in testparm else None)
+            if 'speciescenter' in testparm else None,
+            split=testparm['split'] if 'split' in testparm else 1,
+            )
 
     def test_reacnetgen(self, reacnetgen):
         """Test main process of ReacNetGen."""