diff --git a/Simulation/Aivika/Dynamics/Interpolate.hs b/Simulation/Aivika/Dynamics/Interpolate.hs
--- a/Simulation/Aivika/Dynamics/Interpolate.hs
+++ b/Simulation/Aivika/Dynamics/Interpolate.hs
@@ -12,11 +12,22 @@
 --
 
 module Simulation.Aivika.Dynamics.Interpolate
-       (discreteDynamics,
+       (initDynamics,
+        discreteDynamics,
         interpolateDynamics) where
 
 import Simulation.Aivika.Internal.Specs
 import Simulation.Aivika.Internal.Dynamics
+
+-- | Return the initial value.
+initDynamics :: Dynamics a -> Dynamics a
+{-# INLINE initDynamics #-}
+initDynamics (Dynamics m) =
+  Dynamics $ \p ->
+  let sc = pointSpecs p
+  in m $ p { pointTime = basicTime sc 0 0,
+             pointIteration = 0,
+             pointPhase = 0 }
 
 -- | Discretize the computation in the integration time points.
 discreteDynamics :: Dynamics a -> Dynamics a
diff --git a/Simulation/Aivika/Internal/Parameter.hs b/Simulation/Aivika/Internal/Parameter.hs
--- a/Simulation/Aivika/Internal/Parameter.hs
+++ b/Simulation/Aivika/Internal/Parameter.hs
@@ -48,7 +48,7 @@
 import Control.Monad.Fix
 
 import Data.IORef
-import qualified Data.Map as M
+import qualified Data.IntMap as M
 import Data.Array
 
 import Simulation.Aivika.Generator
diff --git a/Simulation/Aivika/SystemDynamics.hs b/Simulation/Aivika/SystemDynamics.hs
--- a/Simulation/Aivika/SystemDynamics.hs
+++ b/Simulation/Aivika/SystemDynamics.hs
@@ -119,13 +119,13 @@
 --
 
 integEuler :: Dynamics Double
-             -> Double 
              -> Dynamics Double 
+             -> Dynamics Double 
              -> Point -> IO Double
-integEuler (Dynamics f) i (Dynamics y) p = 
+integEuler (Dynamics f) (Dynamics i) (Dynamics y) p = 
   case pointIteration p of
-    0 ->
-      return i
+    0 -> 
+      i p
     n -> do 
       let sc = pointSpecs p
           ty = basicTime sc (n - 1) 0
@@ -136,14 +136,14 @@
       return v
 
 integRK2 :: Dynamics Double
-           -> Double
            -> Dynamics Double
+           -> Dynamics Double
            -> Point -> IO Double
-integRK2 (Dynamics f) i (Dynamics y) p =
+integRK2 (Dynamics f) (Dynamics i) (Dynamics y) p =
   case pointPhase p of
     0 -> case pointIteration p of
       0 ->
-        return i
+        i p
       n -> do
         let sc = pointSpecs p
             ty = basicTime sc (n - 1) 0
@@ -172,14 +172,14 @@
       error "Incorrect phase: integRK2"
 
 integRK4 :: Dynamics Double
-           -> Double
            -> Dynamics Double
+           -> Dynamics Double
            -> Point -> IO Double
-integRK4 (Dynamics f) i (Dynamics y) p =
+integRK4 (Dynamics f) (Dynamics i) (Dynamics y) p =
   case pointPhase p of
     0 -> case pointIteration p of
       0 -> 
-        return i
+        i p
       n -> do
         let sc = pointSpecs p
             ty = basicTime sc (n - 1) 0
@@ -251,11 +251,8 @@
 --           kb = 1
 --       runDynamicsInStopTime $ sequence [a, b, c]
 -- @
---
--- To receive the initial value for an abitrary 'Dynamics' computation,
--- you can always use the 'runDynamicsInStartTime' function.
 integ :: Dynamics Double                  -- ^ the derivative
-         -> Double                        -- ^ the initial value
+         -> Dynamics Double               -- ^ the initial value
          -> Simulation (Dynamics Double)  -- ^ the integral
 integ diff i =
   mdo y <- MU.memoDynamics z
@@ -278,7 +275,7 @@
 -- @     
 smoothI :: Dynamics Double                  -- ^ the value to smooth over time
            -> Dynamics Double               -- ^ time
-           -> Double                        -- ^ the initial value
+           -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth
 smoothI x t i =
   mdo y <- integ ((x - y) / t) i
@@ -291,9 +288,7 @@
 smooth :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth
-smooth x t =
-  do i <- runDynamicsInStartTime x
-     smoothI x t i
+smooth x t = smoothI x t x
 
 -- | Return the third order exponential smooth.
 --
@@ -310,7 +305,7 @@
 -- @     
 smooth3I :: Dynamics Double                  -- ^ the value to smooth over time
             -> Dynamics Double               -- ^ time
-            -> Double                        -- ^ the initial value
+            -> Dynamics Double               -- ^ the initial value
             -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth
 smooth3I x t i =
   mdo y  <- integ ((s2 - y) / t') i
@@ -326,9 +321,7 @@
 smooth3 :: Dynamics Double                  -- ^ the value to smooth over time
            -> Dynamics Double               -- ^ time
            -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth
-smooth3 x t =
-  do i <- runDynamicsInStartTime x
-     smooth3I x t i
+smooth3 x t = smooth3I x t x
 
 -- | Return the n'th order exponential smooth.
 --
@@ -339,7 +332,7 @@
 smoothNI :: Dynamics Double                  -- ^ the value to smooth over time
             -> Dynamics Double               -- ^ time
             -> Int                           -- ^ the order
-            -> Double                        -- ^ the initial value
+            -> Dynamics Double               -- ^ the initial value
             -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth
 smoothNI x t n i =
   mdo s <- forM [1 .. n] $ \k ->
@@ -358,9 +351,7 @@
            -> Dynamics Double               -- ^ time
            -> Int                           -- ^ the order
            -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth
-smoothN x t n =
-  do i <- runDynamicsInStartTime x
-     smoothNI x t n i
+smoothN x t n = smoothNI x t n x
 
 -- | Return the first order exponential delay.
 --
@@ -369,17 +360,15 @@
 --
 -- @
 -- delay1I x t i =
---   mdo t0 <- runDynamicsInStartTime t
---       y  <- integ (x - y \/ t) (i * t0)
+--   mdo y <- integ (x - y \/ t) (i * t)
 --       return $ y \/ t
 -- @     
 delay1I :: Dynamics Double                  -- ^ the value to conserve
            -> Dynamics Double               -- ^ time
-           -> Double                        -- ^ the initial value
+           -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the first order exponential delay
 delay1I x t i =
-  mdo t0 <- runDynamicsInStartTime t
-      y  <- integ (x - y / t) (i * t0)
+  mdo y <- integ (x - y / t) (i * t)
       return $ y / t
 
 -- | Return the first order exponential delay.
@@ -389,20 +378,17 @@
 delay1 :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Simulation (Dynamics Double)  -- ^ the first order exponential delay
-delay1 x t =
-  do i <- runDynamicsInStartTime x
-     delay1I x t i
+delay1 x t = delay1I x t x
 
 -- | Return the third order exponential delay.
 delay3I :: Dynamics Double                  -- ^ the value to conserve
            -> Dynamics Double               -- ^ time
-           -> Double                        -- ^ the initial value
+           -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the third order exponential delay
 delay3I x t i =
-  mdo t0' <- runDynamicsInStartTime t'
-      y   <- integ (s2 / t' - y / t') (i * t0')
-      s2  <- integ (s1 / t' - s2 / t') (i * t0')
-      s1  <- integ (x - s1 / t') (i * t0')
+  mdo y  <- integ (s2 / t' - y / t') (i * t')
+      s2 <- integ (s1 / t' - s2 / t') (i * t')
+      s1 <- integ (x - s1 / t') (i * t')
       let t' = t / 3.0
       return $ y / t'         
 
@@ -413,22 +399,19 @@
 delay3 :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Simulation (Dynamics Double)  -- ^ the third order exponential delay
-delay3 x t =
-  do i <- runDynamicsInStartTime x
-     delay3I x t i
+delay3 x t = delay3I x t x
 
 -- | Return the n'th order exponential delay.
 delayNI :: Dynamics Double                  -- ^ the value to conserve
            -> Dynamics Double               -- ^ time
            -> Int                           -- ^ the order
-           -> Double                        -- ^ the initial value
+           -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay
 delayNI x t n i =
-  mdo t0' <- runDynamicsInStartTime t'
-      s   <- forM [1 .. n] $ \k ->
+  mdo s <- forM [1 .. n] $ \k ->
         if k == 1
-        then integ (x - (a ! 1) / t') (i * t0')
-        else integ ((a ! (k - 1)) / t' - (a ! k) / t') (i * t0')
+        then integ (x - (a ! 1) / t') (i * t')
+        else integ ((a ! (k - 1)) / t' - (a ! k) / t') (i * t')
       let a  = listArray (1, n) s
           t' = t / fromIntegral n
       return $ (a ! n) / t'
@@ -441,9 +424,7 @@
           -> Dynamics Double               -- ^ time
           -> Int                           -- ^ the order
           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay
-delayN x t n =
-  do i <- runDynamicsInStartTime x
-     delayNI x t n i
+delayN x t n = delayNI x t n x
 
 -- | Return the forecast.
 --
@@ -468,20 +449,16 @@
 --
 -- @
 -- trend x at i =
---   mdo x0  <- runDynamicsInStartTime x
---       at0 <- runDynamicsInStartTime at
---       y   <- smoothI x at (x0 \/ (1.0 + i * at0))
---       return $ (x \/ y - 1.0) \/ at
+--   do y <- smoothI x at (x \/ (1.0 + i * at))
+--      return $ (x \/ y - 1.0) \/ at
 -- @
 trend :: Dynamics Double                  -- ^ the value for which the trend is calculated
          -> Dynamics Double               -- ^ the average time
-         -> Double                        -- ^ the initial value
+         -> Dynamics Double               -- ^ the initial value
          -> Simulation (Dynamics Double)  -- ^ the fractional change rate
 trend x at i =
-  mdo x0  <- runDynamicsInStartTime x
-      at0 <- runDynamicsInStartTime at
-      y   <- smoothI x at (x0 / (1.0 + i * at0))
-      return $ (x / y - 1.0) / at
+  do y <- smoothI x at (x / (1.0 + i * at))
+     return $ (x / y - 1.0) / at
 
 --
 -- Difference Equations
@@ -494,15 +471,14 @@
 -- As usual, to create a loopback, you should use the recursive do-notation.
 diffsum :: (Num a, Unboxed a)
            => Dynamics a               -- ^ the difference
-           -> a                        -- ^ the initial value
+           -> Dynamics a               -- ^ the initial value
            -> Simulation (Dynamics a)  -- ^ the sum
-diffsum (Dynamics diff) i =
+diffsum (Dynamics diff) (Dynamics i) =
   mdo y <-
         MU.memo0Dynamics $
         Dynamics $ \p ->
         case pointIteration p of
-          0 ->
-            return i
+          0 -> i p
           n -> do 
             let Dynamics m = y
                 sc = pointSpecs p
@@ -569,9 +545,9 @@
 -- Because of the latter, it allows creating a loop back.
 delayI :: Dynamics a          -- ^ the value to delay
           -> Dynamics Double  -- ^ the lag time
-          -> a                -- ^ the initial value
+          -> Dynamics a       -- ^ the initial value
           -> Simulation (Dynamics a)    -- ^ the delayed value
-delayI (Dynamics x) (Dynamics d) i = M.memo0Dynamics $ Dynamics r 
+delayI (Dynamics x) (Dynamics d) (Dynamics i) = M.memo0Dynamics $ Dynamics r 
   where
     r p = do 
       let t  = pointTime p
@@ -580,7 +556,9 @@
       a <- d p
       let t' = t - a
           n' = fromIntegral $ floor $ (t' - spcStartTime sc) / spcDT sc
-          y | n' < 0    = return i
+          y | n' < 0    = i $ p { pointTime = spcStartTime sc,
+                                  pointIteration = 0, 
+                                  pointPhase = 0 }
             | n' < n    = x $ p { pointTime = t',
                                   pointIteration = n',
                                   pointPhase = -1 }
@@ -604,18 +582,18 @@
 -- @
 -- npv stream rate init factor =
 --   mdo let dt' = liftParameter dt
---       df    <- integ (- df * rate) 1
+--       df <- integ (- df * rate) 1
 --       accum <- integ (stream * df) init
 --       return $ (accum + dt' * stream * df) * factor
 -- @
 npv :: Dynamics Double                  -- ^ the stream
        -> Dynamics Double               -- ^ the discount rate
-       -> Double                        -- ^ the initial value
+       -> Dynamics Double               -- ^ the initial value
        -> Dynamics Double               -- ^ factor
        -> Simulation (Dynamics Double)  -- ^ the Net Present Value (NPV)
 npv stream rate init factor =
   mdo let dt' = liftParameter dt
-      df    <- integ (- df * rate) 1
+      df <- integ (- df * rate) 1
       accum <- integ (stream * df) init
       return $ (accum + dt' * stream * df) * factor
 
@@ -627,22 +605,18 @@
 -- @
 -- npve stream rate init factor =
 --   mdo let dt' = liftParameter dt
---       rate0 <- runDynamicsInStartTime rate
---       dt0   <- liftParameter dt
---       df    <- integ (- df * rate \/ (1 + rate * dt')) (1 \/ (1 + rate0 * dt0))
+--       df <- integ (- df * rate \/ (1 + rate * dt')) (1 \/ (1 + rate * dt'))
 --       accum <- integ (stream * df) init
 --       return $ (accum + dt' * stream * df) * factor
 -- @
 npve :: Dynamics Double                  -- ^ the stream
         -> Dynamics Double               -- ^ the discount rate
-        -> Double                        -- ^ the initial value
+        -> Dynamics Double               -- ^ the initial value
         -> Dynamics Double               -- ^ factor
         -> Simulation (Dynamics Double)  -- ^ the Net Present Value End (NPVE)
 npve stream rate init factor =
   mdo let dt' = liftParameter dt
-      rate0 <- runDynamicsInStartTime rate
-      dt0   <- liftParameter dt
-      df    <- integ (- df * rate / (1 + rate * dt')) (1 / (1 + rate0 * dt0))
+      df <- integ (- df * rate / (1 + rate * dt')) (1 / (1 + rate * dt'))
       accum <- integ (stream * df) init
       return $ (accum + dt' * stream * df) * factor
 
diff --git a/aivika.cabal b/aivika.cabal
--- a/aivika.cabal
+++ b/aivika.cabal
@@ -1,5 +1,5 @@
 name:            aivika
-version:         1.2
+version:         1.2.1
 synopsis:        A multi-paradigm simulation library
 description:
     Aivika is a multi-paradigm simulation library with a strong emphasis
@@ -32,6 +32,8 @@
       operate on the infinite streams of data, because of which some models
       can remind of their high-level graphical representation on the
       diagram used by visual simulation software tools;
+    .
+    * allows simulating circuits with recursive links and delays;
     .
     * supports the activity-oriented paradigm of DES;
     .