diff --git a/.travis_build_linux.sh b/.travis_build_linux.sh
index a9397370877b9b6c252b916035575925aad23b25..4906319cd30c0c9379a32974078b175b4c69a798 100755
--- a/.travis_build_linux.sh
+++ b/.travis_build_linux.sh
@@ -4,3 +4,11 @@ set -e -x
 PATH=/usr/bin:/usr/local/bin:$PATH
 gbp buildpackage  --git-ignore-branch --git-ignore-new -uc -us -d | tee _gbp.log
 pwd 
+ls ../*.deb
+
+if [ -n "$TRAVIS" ]; then
+    echo "We are on travis"
+    sudo dpkg -i ../*.deb
+    sudo apt-get install -f 
+    python -c 'import moose; print(moose.__file__); print(moose.__version__)'
+fi
diff --git a/.travis_build_osx.sh b/.travis_build_osx.sh
index 3f08edd1c586f72efd1b3264c8a1d8ba033c7d01..10d39acf0918246a4b7c185436105355c96b9a5b 100755
--- a/.travis_build_osx.sh
+++ b/.travis_build_osx.sh
@@ -3,4 +3,3 @@ brew tap BhallaLab/moose
 brew install --HEAD moose-nightly
 python -m pip install networkx python-libsbml pyNeuroML  matplotlib
 python -c 'import moose; print( moose.__version__ )'
-python -c 'import moose; print( moose.test() )'
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 756c329d636ea3d03d81c8802eae2d46c89aaa72..cf83ee350e292d81f3dee2513c0b858f9f40d5c2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -108,10 +108,7 @@ if(WITH_GUI)
     add_dependencies(examples  gui)
 endif()
 
-
-# NOTE: leading / is important other whole path will be installed. Do not
-# install usr since it comes with CMAKE_INSTALL_PREFIX.
-install(DIRECTORY ${PYMOOSE_INSTALL_DIR}/usr/
+install(DIRECTORY ${PYMOOSE_INSTALL_DIR}/
     DESTINATION ${CMAKE_INSTALL_PREFIX}
     CONFIGURATIONS Release
     PATTERN ".git" EXCLUDE
diff --git a/moose-core/.travis.yml b/moose-core/.travis.yml
index a89ee791337440c7828cccb68c5424a341c5a962..01188d6d9b71ef4bd491190ad1d7ac71c84e0d90 100644
--- a/moose-core/.travis.yml
+++ b/moose-core/.travis.yml
@@ -1,7 +1,7 @@
 language: cpp
 sudo: required
-group: edge
-
+group: travis_lts
+dist: xenial
 
 os:
 - linux
diff --git a/moose-core/.travis/travis_build_linux.sh b/moose-core/.travis/travis_build_linux.sh
index f715e6c4e1209d0d4f6883d01e0fa449c65afa3d..6885c7c61dc43cbf34ac4d07bae34ccd9bfca3b4 100755
--- a/moose-core/.travis/travis_build_linux.sh
+++ b/moose-core/.travis/travis_build_linux.sh
@@ -40,24 +40,34 @@ unset PYTHONPATH
 $PYTHON2 -m compileall -q .
 if type $PYTHON3 > /dev/null; then $PYTHON3 -m compileall -q . ; fi
 
-# This is only applicable on linux build.
+# Python2 and GSL.
+echo "Currently in `pwd`"
+(
+    mkdir -p _GSL_BUILD && cd _GSL_BUILD
+    cmake -DDEBUG=ON -DPYTHON_EXECUTABLE="$PYTHON2" ..
+    $MAKE && ctest --output-on-failure
+    sudo make install && cd  /tmp
+    $PYTHON2 -c 'import moose;print(moose.__file__);print(moose.version())'
+)
+
+# Python3 with GSL and BOOST. This should be the only build after we drop
+# python2 support.
 echo "Python3: Removed python2-networkx and install python3"
 if type $PYTHON3 > /dev/null; then
     sudo apt-get remove -qq python-networkx || echo "Error with apt"
     sudo apt-get install -qq python3-networkx || echo "Error with apt"
-
     # GSL.
     (
         mkdir -p _GSL_BUILD2 && cd _GSL_BUILD2 && \
             cmake -DPYTHON_EXECUTABLE="$PYTHON3" -DDEBUG=ON ..
-        $MAKE && ctest --output-on-failure 
+        $MAKE && ctest --output-on-failure
     )
 
     # BOOST
     (
         mkdir -p _BOOST_BUILD2 && cd _BOOST_BUILD2 && \
             cmake -DWITH_BOOST_ODE=ON -DPYTHON_EXECUTABLE="$PYTHON3" ..
-        $MAKE && ctest --output-on-failure 
+        $MAKE && ctest --output-on-failure
     )
     echo "All done"
 else
diff --git a/moose-core/.travis/travis_prepare_linux.sh b/moose-core/.travis/travis_prepare_linux.sh
index 484f1a4ec09e5b4ac3f4a4e3b778a5529746d00f..0c0e604f3588ba294f8c86bcdf99142ca9e17884 100755
--- a/moose-core/.travis/travis_prepare_linux.sh
+++ b/moose-core/.travis/travis_prepare_linux.sh
@@ -10,7 +10,7 @@
 #       OPTIONS: ---
 #  REQUIREMENTS: ---
 #          BUGS: ---
-#         NOTES: ---
+#         NOTES: Always run with sudo permission.
 #        AUTHOR: Dilawar Singh (), dilawars@ncbs.res.in
 #  ORGANIZATION: NCBS Bangalore
 #       CREATED: 01/02/2017 10:10:02 AM
@@ -21,15 +21,24 @@ set -o nounset                              # Treat unset variables as an error
 set +e  # Let installation fail in some command
 
 apt-get install -qq libxml2-dev libbz2-dev
-apt-get install -qq libhdf5-serial-dev
 apt-get install -qq make cmake
 apt-get install -qq python-numpy python-matplotlib python-networkx python-pip
+apt-get install -qq python3-lxml python-lxml
 apt-get install -qq python3-numpy python3-matplotlib python3-dev
-apt-get install -qq libboost-all-dev
-apt-get install -qq libgsl0-dev
 apt-get install -qq python-pip python3-pip
 apt-get install -qq libgraphviz-dev
 
+# Gsl
+apt-get install -qq libgsl0-dev || apt-get install -qq libgsl-dev
+
+# Boost related.
+apt-get install -qq liblapack-dev
+apt-get install -qq libboost-all-dev
+
 # Dependencies for NML2
 apt-get install -qq python-scipy python3-scipy
-#pip install pyNeuroML libNeuroML
+apt-get install -qq python-lxml python3-lxml
+apt-get install -qq python-setuptools python3-setuptools
+apt-get install -qq python-tornado python3-tornado
+python2 -m pip install pyNeuroML libNeuroML
+python3 -m pip install pyNeuroML libNeuroML
diff --git a/moose-core/CMakeLists.txt b/moose-core/CMakeLists.txt
index 3be4235b3d563ff57adfe0d952ebecc84c44d993..a71896f639068f44b6c42b11c73be926c7253a80 100644
--- a/moose-core/CMakeLists.txt
+++ b/moose-core/CMakeLists.txt
@@ -76,6 +76,8 @@ else()
 endif()
 
 ################################ CMAKE OPTIONS ##################################
+option(WITH_NSDF            "Enable NSDF support. Requires hdf5" OFF )
+
 option(DEBUG                "Build with debug support" OFF)
 option(GPROF                "Build for profiling using gprof" OFF)
 option(ENABLE_UNIT_TESTS    "Enable unit tests (DEBUG should also be ON)" OFF)
@@ -187,68 +189,71 @@ if(WITH_BOOST_ODE)
     add_definitions(-DUSE_BOOST_ODE -UUSE_GSL)
 endif()
 
-
-find_package(HDF5 COMPONENTS CXX HL)
-if(NOT HDF5_FOUND)
-    message(
-        "==================================================================\n"
-        " HDF5 not found. Disabling support. Required for nsdf. \n\n"
-        " If you need hdf5 support, please install hdf5-dev or hdf5-devel\n"
-        " package or equivalent.\n\n"
-        "     $ sudo apt-get install libhdf5-dev \n"
-        "     $ sudo yum install libhdf5-devel \n"
-        "     $ brew install hdf5 \n\n"
-        " Otherwise, continue with 'make' and 'make install' \n"
-        " If you install hdf5 to non-standard path, export environment \n"
-        " variable HDF5_ROOT to the location. Rerun cmake \n"
-        "================================================================ \n"
-        )
-endif(NOT HDF5_FOUND)
-
-if(HDF5_FOUND)
-    include_directories( ${HDF5_INCLUDE_DIRS} )
-    add_definitions( -DUSE_HDF5 )
-    if(HDF5_USE_STATIC_LIBRARIES)
-	message(STATUS "Finding static HDF5 libraries in $ENV{HDF5_ROOT}")
-        find_library(HDF5_CXX_LIBRARIES NAMES libhdf5.a
-            PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
+# NSDF5 support. Disabled by default.
+if(WITH_NSDF)
+    find_package(HDF5 COMPONENTS CXX HL)
+    if(NOT HDF5_FOUND)
+        message(
+            "==================================================================\n"
+            " HDF5 not found. Disabling NSDF support.\n\n"
+            " If you need NSDF support, please install hdf5-dev or hdf5-devel\n"
+            " package or equivalent.\n\n"
+            "     $ sudo apt-get install libhdf5-dev \n"
+            "     $ sudo yum install libhdf5-devel \n"
+            "     $ brew install hdf5 \n\n"
+            " Otherwise, continue with 'make' and 'make install' \n"
+            " If you install hdf5 to non-standard path, export environment \n"
+            " variable HDF5_ROOT to the location. Rerun cmake \n"
+            "================================================================ \n"
             )
-        find_library(HDF5_HL_LIBRARIES NAMES libhdf5_hl.a
-            PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
-            )
-        set(HDF5_LIBRARIES ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES})
+    endif(NOT HDF5_FOUND)
+
+    if(HDF5_FOUND)
+        include_directories( ${HDF5_INCLUDE_DIRS} )
+        add_definitions( -DUSE_HDF5 -DENABLE_NSDF )
+        if(HDF5_USE_STATIC_LIBRARIES)
+            message(STATUS "Finding static HDF5 libraries in $ENV{HDF5_ROOT}")
+            find_library(HDF5_CXX_LIBRARIES NAMES libhdf5.a
+                PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
+                )
+            find_library(HDF5_HL_LIBRARIES NAMES libhdf5_hl.a
+                PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
+                )
+            set(HDF5_LIBRARIES ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES})
+        endif()
+
+
+        # Make sure, HDF5_HL_LIBRARIES are set. The COMPONENTS in find_package may
+        # or may not work. See BhallaLab/moose-core#163.
+        if(NOT HDF5_HL_LIBRARIES)
+            set(HDF5_HL_LIBRARIES ${HDF5_HL_LIBRARIES})
+        endif(NOT HDF5_HL_LIBRARIES)
+        list(APPEND HDF5_LIBRARIES ${HDF5_HL_LIBRARIES})
+
+        # message(STATUS "MOOSE will use following HDF5 ${HDF5_LIBRARIES}" )
+        foreach(HDF5_LIB ${HDF5_LIBRARIES})
+            if(HDF5_LIB)
+                get_filename_component( HDF5_LIB_EXT ${HDF5_LIB} EXT )
+                if(HDF5_LIB_EXT)
+                    if(${HDF5_LIB_EXT} STREQUAL ".a")
+                        list(APPEND STATIC_LIBRARIES ${HDF5_LIB} )
+                    else( )
+                        list(APPEND SYSTEM_SHARED_LIBS ${HDF5_LIB} )
+                    endif( )
+                endif()
+            endif( )
+        endforeach( )
+    else( HDF5_FOUND )
+        message(STATUS "HDF5 is not found" )
+    endif( HDF5_FOUND )
+    # This is a fix for new HDF5 package on Debian/Ubuntu which installs hdf5
+    # headers in non-standard path. issue #80.
+    if(HDF5_LIBRARY_DIRS)
+        set_target_properties( libmoose PROPERTIES LINK_FLAGS "-L${HDF5_LIBRARY_DIRS}" )
     endif()
-
-
-    # Make sure, HDF5_HL_LIBRARIES are set. The COMPONENTS in find_package may
-    # or may not work. See BhallaLab/moose-core#163.
-    if(NOT HDF5_HL_LIBRARIES)
-        set(HDF5_HL_LIBRARIES ${HDF5_HL_LIBRARIES})
-    endif(NOT HDF5_HL_LIBRARIES)
-    list(APPEND HDF5_LIBRARIES ${HDF5_HL_LIBRARIES})
-
-    # message(STATUS "MOOSE will use following HDF5 ${HDF5_LIBRARIES}" )
-    foreach(HDF5_LIB ${HDF5_LIBRARIES})
-        if(HDF5_LIB)
-            get_filename_component( HDF5_LIB_EXT ${HDF5_LIB} EXT )
-            if(HDF5_LIB_EXT)
-                if(${HDF5_LIB_EXT} STREQUAL ".a")
-                    list(APPEND STATIC_LIBRARIES ${HDF5_LIB} )
-                else( )
-                    list(APPEND SYSTEM_SHARED_LIBS ${HDF5_LIB} )
-                endif( )
-            endif()
-        endif( )
-    endforeach( )
-else( HDF5_FOUND )
-    message(STATUS "HDF5 is not found" )
-endif( HDF5_FOUND )
-
-# This is a fix for new HDF5 package on Debian/Ubuntu which installs hdf5
-# headers in non-standard path. issue #80.
-if(HDF5_LIBRARY_DIRS)
-    set_target_properties( libmoose PROPERTIES LINK_FLAGS "-L${HDF5_LIBRARY_DIRS}" )
-endif()
+else(WITH_NSDF)
+    message(STATUS "NSDF support is disabled" )
+endif(WITH_NSDF)
 
 # Openmpi
 if(WITH_MPI)
@@ -432,7 +437,7 @@ message(STATUS "binary distribution file ${PYMOOSE_BDIST_FILE}")
 add_custom_target(bdist ALL DEPENDS  ${PYMOOSE_BDIST_FILE} )
 
 # Any command using setup.cmake.py must run in the same directory. Use option
-# `--relative` to prefix is aways fixed. 
+# `--relative` to prefix is aways fixed.
 add_custom_command( OUTPUT ${PYMOOSE_BDIST_FILE}
     COMMAND ${PYTHON_EXECUTABLE} setup.cmake.py bdist_dumb -p ${_platform}
         -d ${PYMOOSE_BDIST_DIR} --relative
@@ -464,7 +469,7 @@ add_dependencies( copy_pymoose _moose )
 install(TARGETS moose.bin DESTINATION bin CONFIGURATIONS Debug)
 install(TARGETS libmoose DESTINATION lib CONFIGURATIONS Debug)
 
-# install pymoose bdist. 
+# install pymoose bdist.
 install(DIRECTORY ${PYMOOSE_BDIST_INSTALL_DIR}/
     DESTINATION ${CMAKE_INSTALL_PREFIX}
     CONFIGURATIONS Release Debug
diff --git a/moose-core/biophysics/ChanCommon.cpp b/moose-core/biophysics/ChanCommon.cpp
index 9a015e4fc90f38cd64b9b9dc5e42840f68c31c3b..124fe67f7e5bb337cc3dbf8d15b2f6c786f365e5 100644
--- a/moose-core/biophysics/ChanCommon.cpp
+++ b/moose-core/biophysics/ChanCommon.cpp
@@ -15,17 +15,19 @@
 // Constructor
 ///////////////////////////////////////////////////
 ChanCommon::ChanCommon()
-			:
-			Vm_( 0.0 ),
-			Gbar_( 0.0 ), modulation_( 1.0 ),
-			Ek_( 0.0 ),
-			Gk_( 0.0 ), Ik_( 0.0 )
+    :
+    Vm_( 0.0 ),
+    Gbar_( 0.0 ), modulation_( 1.0 ),
+    Ek_( 0.0 ),
+    Gk_( 0.0 ), Ik_( 0.0 )
 {
-	;
+    ;
 }
 
 ChanCommon::~ChanCommon()
-{;}
+{
+    ;
+}
 
 ///////////////////////////////////////////////////
 // Field function definitions
@@ -33,54 +35,54 @@ ChanCommon::~ChanCommon()
 
 void ChanCommon::vSetGbar( const Eref& e, double Gbar )
 {
-	Gbar_ = Gbar;
+    Gbar_ = Gbar;
 }
 
 double ChanCommon::vGetGbar( const Eref& e ) const
 {
-	return Gbar_;
+    return Gbar_;
 }
 
 void ChanCommon::vSetModulation( const Eref& e, double modulation )
 {
-	modulation_ = modulation;
+    modulation_ = modulation;
 }
 
 double ChanCommon::vGetModulation( const Eref& e ) const
 {
-	return modulation_;
+    return modulation_;
 }
 
 double ChanCommon::getModulation() const
 {
-	return modulation_;
+    return modulation_;
 }
 
 void ChanCommon::vSetEk( const Eref& e, double Ek )
 {
-	Ek_ = Ek;
+    Ek_ = Ek;
 }
 double ChanCommon::vGetEk( const Eref& e ) const
 {
-	return Ek_;
+    return Ek_;
 }
 
 void ChanCommon::vSetGk( const Eref& e, double Gk )
 {
-	Gk_ = Gk;
+    Gk_ = Gk;
 }
 double ChanCommon::vGetGk( const Eref& e ) const
 {
-	return Gk_;
+    return Gk_;
 }
 
 void ChanCommon::vSetIk( const Eref& e, double Ik )
 {
-	Ik_ = Ik;
+    Ik_ = Ik;
 }
 double ChanCommon::vGetIk( const Eref& e ) const
 {
-	return Ik_;
+    return Ik_;
 }
 
 ///////////////////////////////////////////////////
@@ -89,7 +91,7 @@ double ChanCommon::vGetIk( const Eref& e ) const
 
 void ChanCommon::vHandleVm( double Vm )
 {
-	Vm_ = Vm;
+    Vm_ = Vm;
 }
 
 ///////////////////////////////////////////////////
@@ -99,34 +101,34 @@ void ChanCommon::vHandleVm( double Vm )
 
 void ChanCommon::sendProcessMsgs(  const Eref& e, const ProcPtr info )
 {
-		ChanBase::channelOut()->send( e, Gk_, Ek_ );
-	// This is used if the channel connects up to a conc pool and
-	// handles influx of ions giving rise to a concentration change.
-		ChanBase::IkOut()->send( e, Ik_ );
-	// Needed by GHK-type objects
-		ChanBase::permeability()->send( e, Gk_ );
+    ChanBase::channelOut()->send( e, Gk_, Ek_ );
+    // This is used if the channel connects up to a conc pool and
+    // handles influx of ions giving rise to a concentration change.
+    ChanBase::IkOut()->send( e, Ik_ );
+    // Needed by GHK-type objects
+    ChanBase::permeability()->send( e, Gk_ );
 }
 
 
 void ChanCommon::sendReinitMsgs(  const Eref& e, const ProcPtr info )
 {
-		ChanBase::channelOut()->send( e, Gk_, Ek_ );
-		ChanBase::IkOut()->send( e, Ik_ );
-	// Needed by GHK-type objects
-		ChanBase::permeability()->send( e, Gk_ );
+    ChanBase::channelOut()->send( e, Gk_, Ek_ );
+    ChanBase::IkOut()->send( e, Ik_ );
+    // Needed by GHK-type objects
+    ChanBase::permeability()->send( e, Gk_ );
 }
 
 void ChanCommon::updateIk()
 {
-	Ik_ = ( Ek_ - Vm_ ) * Gk_;
+    Ik_ = ( Ek_ - Vm_ ) * Gk_;
 }
 
 double ChanCommon::getVm() const
 {
-	return Vm_;
+    return Vm_;
 }
 
 double ChanCommon::getGbar() const
 {
-	return Gbar_;
+    return Gbar_;
 }
diff --git a/moose-core/biophysics/ChanCommon.h b/moose-core/biophysics/ChanCommon.h
index 007c2f3785197cbb942224a1a475213c8dc0b536..50725bb6679114534852e0e010612b180f1613b3 100644
--- a/moose-core/biophysics/ChanCommon.h
+++ b/moose-core/biophysics/ChanCommon.h
@@ -20,78 +20,75 @@
 
 class ChanCommon: public virtual ChanBase
 {
-	public:
-		ChanCommon();
-		~ChanCommon();
-
-		/////////////////////////////////////////////////////////////
-		// Value field access function definitions
-		/////////////////////////////////////////////////////////////
-
-		void vSetGbar( const Eref& e, double Gbar );
-		double vGetGbar( const Eref& e ) const;
-		void vSetModulation( const Eref& e, double modulation );
-		double vGetModulation( const Eref& e ) const;
-		double getModulation() const;
-		void vSetEk( const Eref& e, double Ek );
-		double vGetEk( const Eref& e ) const;
-		void vSetGk( const Eref& e, double Gk );
-		double vGetGk( const Eref& e ) const;
-		/// Ik is read-only for MOOSE, but we provide the set
-		/// func for derived classes to update it.
-		void vSetIk( const Eref& e, double Ic );
-		double vGetIk( const Eref& e ) const;
-
-		/////////////////////////////////////////////////////////////
-		// Dest function definitions
-		/////////////////////////////////////////////////////////////
-
-		/**
-		 * Assign the local Vm_ to the incoming Vm from the compartment
-		 */
-		void vHandleVm( double Vm );
-
-		/////////////////////////////////////////////////////////////
-		/**
-		 * This function sends out the messages expected of a channel,
-		 * after process. Used as a utility by various derived classes.
-		 */
-		void sendProcessMsgs( const Eref& e, const ProcPtr info );
-		void sendReinitMsgs( const Eref& e, const ProcPtr info );
-
-		/////////////////////////////////////////////////////////////
-
-		/**
-		 * Utility function for a common computation using local variables
-		 */
-		void updateIk();
-
-
-		/// Utility function to access Vm
-		double getVm() const;
-
-		/// Utility function to acces Gbar
-		double getGbar() const;
-
-		/// Specify the Class Info static variable for initialization.
-		static const Cinfo* initCinfo();
-	protected:
-		/// Vm_ is input variable from compartment, used for most rates
-		double Vm_;
-
-	private:
-		/// Channel maximal conductance
-		double Gbar_;
-		/// Channel modulation. Scales conductance.
-		double modulation_;
-		/// Reversal potential of channel
-		double Ek_;
-
-		/// Channel actual conductance depending on opening of gates.
-		double Gk_;
-		/// Channel current
-		double Ik_;
+public:
+    ChanCommon();
+    ~ChanCommon();
+
+    /////////////////////////////////////////////////////////////
+    // Value field access function definitions
+    /////////////////////////////////////////////////////////////
+
+    void vSetGbar( const Eref& e, double Gbar );
+    double vGetGbar( const Eref& e ) const;
+    void vSetModulation( const Eref& e, double modulation );
+    double vGetModulation( const Eref& e ) const;
+    double getModulation() const;
+    void vSetEk( const Eref& e, double Ek );
+    double vGetEk( const Eref& e ) const;
+    void vSetGk( const Eref& e, double Gk );
+    double vGetGk( const Eref& e ) const;
+    /// Ik is read-only for MOOSE, but we provide the set
+    /// func for derived classes to update it.
+    void vSetIk( const Eref& e, double Ic );
+    double vGetIk( const Eref& e ) const;
+
+    /////////////////////////////////////////////////////////////
+    // Dest function definitions
+    /////////////////////////////////////////////////////////////
+
+    /**
+     * Assign the local Vm_ to the incoming Vm from the compartment
+     */
+    void vHandleVm( double Vm );
+
+    /////////////////////////////////////////////////////////////
+    /**
+     * This function sends out the messages expected of a channel,
+     * after process. Used as a utility by various derived classes.
+     */
+    void sendProcessMsgs( const Eref& e, const ProcPtr info );
+    void sendReinitMsgs( const Eref& e, const ProcPtr info );
+
+    /////////////////////////////////////////////////////////////
+
+    /**
+     * Utility function for a common computation using local variables
+     */
+    void updateIk();
+
+    /// Utility function to access Vm
+    double getVm() const;
+
+    /// Utility function to acces Gbar
+    double getGbar() const;
+
+    /// Specify the Class Info static variable for initialization.
+    static const Cinfo* initCinfo();
+protected:
+    /// Vm_ is input variable from compartment, used for most rates
+    double Vm_;
+
+private:
+    /// Channel maximal conductance
+    double Gbar_;
+    /// Channel modulation. Scales conductance.
+    double modulation_;
+    /// Reversal potential of channel
+    double Ek_;
+
+    /// Channel actual conductance depending on opening of gates.
+    double Gk_;
+    /// Channel current
+    double Ik_;
 };
-
-
 #endif // _ChanCommon_h
diff --git a/moose-core/biophysics/Leakage.cpp b/moose-core/biophysics/Leakage.cpp
index ae78b6ee2e43de1a4ef83d8f1052c9bd9ad27fe5..b3bf9b71a1ee11cb6278c545e67183f784ed1f85 100644
--- a/moose-core/biophysics/Leakage.cpp
+++ b/moose-core/biophysics/Leakage.cpp
@@ -52,7 +52,8 @@
 
 const Cinfo* Leakage::initCinfo()
 {
-    static string doc[] = {
+    static string doc[] =
+    {
         "Name", "Leakage",
         "Author", "Subhasis Ray, 2009, Upi Bhalla 2014 NCBS",
         "Description", "Leakage: Passive leakage channel."
@@ -63,8 +64,8 @@ const Cinfo* Leakage::initCinfo()
     static Cinfo LeakageCinfo(
         "Leakage",
         ChanBase::initCinfo(),
-		0,
-		0,
+        0,
+        0,
         &dinfo,
         doc,
         sizeof( doc ) / sizeof( string ));
@@ -88,22 +89,22 @@ Leakage::~Leakage()
 
 void Leakage::vProcess( const Eref & e, ProcPtr p )
 {
-	ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
-	updateIk();
+    ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
+    updateIk();
     sendProcessMsgs(e, p);
 }
 
 void Leakage::vReinit( const Eref & e, ProcPtr p )
 {
-	ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
-	updateIk();
+    ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
+    updateIk();
     sendReinitMsgs(e, p);
 }
 
 void Leakage::vSetGbar( const Eref& e, double gbar )
 {
-		ChanCommon::vSetGk( e, gbar * this->vGetModulation( e ) );
-		ChanCommon::vSetGbar( e, gbar );
+    ChanCommon::vSetGk( e, gbar * this->vGetModulation( e ) );
+    ChanCommon::vSetGbar( e, gbar );
 }
 
 //
diff --git a/moose-core/biophysics/Leakage.h b/moose-core/biophysics/Leakage.h
index 959b5327a1e566b4337f2a0164d6196670b45fa7..3b20eaf797b15bfbfa784486f4058a5eca784bc8 100644
--- a/moose-core/biophysics/Leakage.h
+++ b/moose-core/biophysics/Leakage.h
@@ -51,7 +51,7 @@
 
 class Leakage: public ChanCommon
 {
-  public:
+public:
     Leakage();
     ~Leakage();
     void vProcess( const Eref & e, ProcPtr p );
diff --git a/moose-core/builtins/CMakeLists.txt b/moose-core/builtins/CMakeLists.txt
index bbc4ab70e5281c3068213f6716adba046bb49093..64caa56166ed0465785cbedccc1ed120619a6ead 100644
--- a/moose-core/builtins/CMakeLists.txt
+++ b/moose-core/builtins/CMakeLists.txt
@@ -3,7 +3,7 @@ include_directories( ${CMAKE_SOURCE_DIR}/basecode )
 include_directories( ${CMAKE_SOURCE_DIR}/external/muparser/include )
 include_directories( ${CMAKE_SOURCE_DIR}/scheduling )
 
-add_library(moose_builtins
+set(SRCS
     Arith.cpp
     Group.cpp
     Mstring.cpp
@@ -20,10 +20,17 @@ add_library(moose_builtins
     Streamer.cpp
     Stats.cpp
     Interpol2D.cpp
-    HDF5WriterBase.cpp
-    NSDFWriter.cpp
-    HDF5DataWriter.cpp
     SpikeStats.cpp
     testBuiltins.cpp
-    testNSDF.cpp
     )
+
+if(WITH_NSDF AND HDF5_FOUND)
+    list(APPEND SRCS
+        HDF5WriterBase.cpp
+        NSDFWriter.cpp
+        HDF5DataWriter.cpp
+        testNSDF.cpp
+        )
+endif()
+
+add_library(moose_builtins ${SRCS} )
diff --git a/moose-core/builtins/testBuiltins.cpp b/moose-core/builtins/testBuiltins.cpp
index c301a6a27331df29739bc53bd600794ce6ac0744..4c05470ae1bcc4e7b859b0ee9d7974890595bb38 100644
--- a/moose-core/builtins/testBuiltins.cpp
+++ b/moose-core/builtins/testBuiltins.cpp
@@ -18,7 +18,9 @@
 
 #include "../shell/Shell.h"
 
+#ifdef ENABLE_NSDF
 extern void testNSDF();
+#endif
 
 void testArith()
 {
@@ -425,7 +427,9 @@ void testBuiltins()
 {
 	testArith();
 	testTable();
+#if ENABLE_NSDF
         testNSDF();
+#endif
 }
 
 void testBuiltinsProcess()
diff --git a/moose-core/ksolve/CMakeLists.txt b/moose-core/ksolve/CMakeLists.txt
index f86c69c08d82ae568b283bacb5f2afde6baa2a93..c43875a011951603ac3e865e9b73ccac7b9f279f 100644
--- a/moose-core/ksolve/CMakeLists.txt
+++ b/moose-core/ksolve/CMakeLists.txt
@@ -2,9 +2,21 @@ cmake_minimum_required(VERSION 2.8)
 include(CheckIncludeFileCXX)
 
 include_directories(../builtins ../basecode ../utility ../kinetics)
-include_directories(../ )
 include_directories(../mesh)
 include_directories(../external/muparser/include )
+
+if(WITH_BOOST)
+    find_package(Boost REQUIRED COMPONENTS thread)
+    add_definitions( -DUSE_BOOST_ASYNC )
+    if("${CMAKE_BUILD_TYPE}" STREQUAL "Release")
+        add_definitions( -UBOOST_UBLAS_TYPE_CHECK )
+    else()
+        add_definitions( -DBOOST_UBLAS_TYPE_CHECK )
+    endif()
+    set(WITH_BOOST_ODE ON)
+    include_directories( ${Boost_INCLUDE_DIRS} )
+endif()
+
 IF(WITH_BOOST_ODE)
     check_include_file_cxx( ${Boost_INCLUDE_DIRS}/boost/numeric/odeint.hpp
         ODEINT_EXISTS)
@@ -18,30 +30,24 @@ elseif(WITH_GSL)
     include_directories( ${GSL_INCLUDE_DIRS} )
 endif(WITH_BOOST_ODE)
 
-if(WITH_BOOST)
-    find_package(Boost REQUIRED COMPONENTS thread)
-    add_definitions( -DUSE_BOOST_ASYNC )
-    include_directories( ${Boost_INCLUDE_DIRS} )
-endif()
-
 set(KSOLVE_SRCS
-	KinSparseMatrix.cpp
-	ZombiePool.cpp
-        ZombieFunction.cpp
-        ZombieBufPool.cpp
-	ZombieReac.cpp
-	ZombieEnz.cpp
-	ZombieMMenz.cpp
-        VoxelPoolsBase.cpp
-	VoxelPools.cpp
-        GssaVoxelPools.cpp
-	RateTerm.cpp
-        FuncTerm.cpp
-	Stoich.cpp
-	Ksolve.cpp
-        Gsolve.cpp
-        ZombiePoolInterface.cpp
-        testKsolve.cpp
+    KinSparseMatrix.cpp
+    ZombiePool.cpp
+    ZombieFunction.cpp
+    ZombieBufPool.cpp
+    ZombieReac.cpp
+    ZombieEnz.cpp
+    ZombieMMenz.cpp
+    VoxelPoolsBase.cpp
+    VoxelPools.cpp
+    GssaVoxelPools.cpp
+    RateTerm.cpp
+    FuncTerm.cpp
+    Stoich.cpp
+    Ksolve.cpp
+    Gsolve.cpp
+    ZombiePoolInterface.cpp
+    testKsolve.cpp
     )
 
 if(WITH_GSL)
diff --git a/moose-core/ksolve/Ksolve.cpp b/moose-core/ksolve/Ksolve.cpp
index 81bc118fc29118f5a746740710879a286468ffea..a072b4c0461110ef058dd7bbe60205df14bef1a2 100644
--- a/moose-core/ksolve/Ksolve.cpp
+++ b/moose-core/ksolve/Ksolve.cpp
@@ -235,11 +235,7 @@ static const Cinfo* ksolveCinfo = Ksolve::initCinfo();
 
 Ksolve::Ksolve()
     :
-#ifdef USE_GSL
     method_( "rk5" ),
-#elif USE_BOOST_ODE
-    method_( "rk5a" ),
-#endif
     epsAbs_( 1e-7 ),
     epsRel_( 1e-7 ),
     numThreads_( 1 ),
diff --git a/moose-core/ksolve/NonlinearSystem.h b/moose-core/ksolve/NonlinearSystem.h
index 5792f28773388b0796aefe980cc065d21d7586b6..e2ea270f645b7480f444fd09521698c87ef06872 100644
--- a/moose-core/ksolve/NonlinearSystem.h
+++ b/moose-core/ksolve/NonlinearSystem.h
@@ -31,8 +31,6 @@
 #include <boost/numeric/ublas/lu.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/io.hpp>
-
-
 #include "VoxelPools.h"
 
 using namespace std;
@@ -46,10 +44,10 @@ typedef ublas::matrix<value_type> matrix_type;
 class ReacInfo
 {
     public:
-        int rank;
-        int num_reacs;
+        size_t rank;
+        size_t num_reacs;
         size_t num_mols;
-        int nIter;
+        size_t nIter;
         double convergenceCriterion;
         double* T;
         VoxelPools* pool;
@@ -60,8 +58,10 @@ class ReacInfo
 
 
 /* Matrix inversion routine.
-   Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
-    template<class T>
+ * Uses lu_factorize and lu_substitute in uBLAS to invert a matrix 
+ */
+
+template<class T>
 bool inverse(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
 {
     using namespace boost::numeric::ublas;
@@ -73,6 +73,7 @@ bool inverse(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
 
     // perform LU-factorization
     int res = lu_factorize(A,pm);
+
     if( res != 0 ) return false;
 
     // create identity matrix of "inverse"
@@ -103,14 +104,19 @@ public:
         x1.resize( size_, 0);
 
         ri.nVec.resize( size_ );
-        dx_ = sqrt( numeric_limits<double>::epsilon() );
+
+        // Find machine epsilon to decide on dx. Machine eps is typically 2.22045e-16 on 64 bit machines/intel. 
+        // Square root of this value is a very good value. The values matches with GSL solver nicely.
+        // Making this value too small or too big gonna cause error in blas methods especially when 
+        // computing inverse.  
+        // Gnu-GSL uses GSL_SQRT_DBL_EPSILON  which is 1.4901161193847656e-08
+        // dx_ = 1.4901161193847656e-08;
+        dx_ = std::sqrt( numeric_limits<double>::epsilon() );
     }
 
-    vector_type compute_at(const vector_type& x)
+    bool compute_at(const vector_type& x, vector_type& result)
     {
-        vector_type result( size_ );
-        system(x, result);
-        return result;
+        return system(x, result);
     }
 
     int apply( )
@@ -118,22 +124,51 @@ public:
         return system(x_, f_);
     }
 
-    int compute_jacobians( const vector_type& x, bool compute_inverse = true )
+    int compute_jacobians( const vector_type& x, bool compute_inverse = false )
     {
         for( size_t i = 0; i < size_; i++)
+        {
+            // This trick I leart by looking at GSL implmentation.
+            double dx = dx_ * std::fabs(x[i]);
+            if( dx == 0 )
+                dx = dx_;
+
             for( size_t j = 0; j < size_; j++)
             {
                 vector_type temp = x;
-                temp[j] += dx_;
-                J_(i, j) = (compute_at(temp)[i] - compute_at(x)[i]) / dx_;
+                temp[j] += dx;
+                vector_type res1, res2;
+                auto r1 = compute_at(temp, res1);
+                auto r2 = compute_at(x, res2);
+
+                if( 0 == r1 && 0 == r2 )
+                {
+                    J_(i, j) = (res1[i] - res2[i]) / dx;
+                    if( std::isnan(J_(i,j)) || std::isinf(J_(i,j)) )
+                    {
+                        /* Try increasing dx */
+                        J_.clear();
+                        return -1;
+                    }
+                }
+                else
+                {
+                    J_.clear();
+                    return -1;
+                }
             }
+        }
 
-        // is_jacobian_valid_ = true;
-        // Keep the inverted J_ ready
-        //if(is_jacobian_valid_ and compute_inverse )
         if( compute_inverse )
-            inverse( J_, invJ_ );
-
+        {
+            try {
+                inverse( J_, invJ_ );
+            } catch ( exception & e ) {
+                J_.clear();
+                invJ_.clear();
+                return -1;
+            }
+        }
         return 0;
     }
 
@@ -147,9 +182,15 @@ public:
             init[i] = x[i];
 
         x_ = init;
-        apply();
-
-        compute_jacobians( init );
+        if( 0 == apply() )
+        {
+            if( 0 != compute_jacobians( init ) )
+            {
+                return;
+            }
+        }
+        else
+            return;
     }
 
     string to_string( )
@@ -167,22 +208,16 @@ public:
 
     int system( const vector_type& x, vector_type& f )
     {
-        int num_consv = ri.num_mols - ri.rank;
+        size_t num_consv = ri.num_mols - ri.rank;
         for ( size_t i = 0; i < ri.num_mols; ++i )
         {
             double temp = x[i] * x[i] ;
 
-#if 0
             // if overflow
-            if ( std::isnan( temp ) or std::isinf( temp ) )
-            {
-                cerr << "Failed: ";
-                for( auto v : ri.nVec ) cerr << v << ", ";
-                cerr << endl;
+            if ( std::isnan( temp ) || std::isinf( temp ) )
                 return -1;
-            }
-#endif
-            ri.nVec[i] = temp;
+            else
+                ri.nVec[i] = temp;
         }
 
         vector< double > vels;
@@ -192,21 +227,21 @@ public:
 
         // y = Nr . v
         // Note that Nr is row-echelon: diagonal and above.
-        for ( int i = 0; i < ri.rank; ++i )
+        f.resize( ri.rank + num_consv);
+        for ( size_t i = 0; i < ri.rank; ++i )
         {
             double temp = 0;
-            for ( int j = i; j < ri.num_reacs; ++j )
+            for ( size_t j = i; j < ri.num_reacs; ++j )
                 temp += ri.Nr(i, j ) * vels[j];
             f[i] = temp ;
         }
 
         // dT = gamma.S - T
-        for ( int i = 0; i < num_consv; ++i )
+        for ( size_t i = 0; i < num_consv; ++i )
         {
             double dT = - ri.T[i];
             for ( size_t  j = 0; j < ri.num_mols; ++j )
                 dT += ri.gamma(i, j) * x[j] * x[j];
-
             f[ i + ri.rank] = dT ;
         }
         return 0;
@@ -222,99 +257,42 @@ public:
      * @return  If successful, return true. Check the variable `x_` at
      * which the system f_ is close to zero (within  the tolerance).
      */
-    bool find_roots_gnewton( double tolerance = 1e-7 , size_t max_iter = 50)
+    bool find_roots_gnewton( double tolerance, size_t max_iter)
     {
-        //tolerance = sqrt( numeric_limits<double>::epsilon() );
         double norm2OfDiff = 1.0;
         size_t iter = 0;
-        int status = apply();
+        if(0 != apply() )
+            return false;
 
-        while( ublas::norm_2(f_) >= tolerance )
+        while( ublas::norm_2(f_) > tolerance )
         {
             iter += 1;
-            compute_jacobians( x_, true );
+            ri.nIter = iter;
+
+            if( 0 != compute_jacobians( x_, true ) )
+            {
+                J_.clear();
+                invJ_.clear();
+                return false;
+            }
+
             vector_type correction = ublas::prod( invJ_, f_ );
             x_ -=  correction;
 
             // If could not compute the value of system successfully.
-            status = apply();
-            if( 0 != status )
+            if( 0 != apply() )
+            {
+                x_.clear();
                 return false;
+            }
 
             if( iter >= max_iter )
                 break;
 
         }
-
-        ri.nIter = iter;
-
-        if( iter >= max_iter )
-            return false;
-
         return true;
     }
 
-    /**
-     * @brief Compute the slope of function in given dimension.
-     *
-     * @param which_dimen The index of dimension.
-     *
-     * @return  Slope.
-     */
-    value_type slope( unsigned int which_dimen )
-    {
-        vector_type x = x_;
-        x[which_dimen] += dx_;
-        // x1 and x2 holds the f_ of system at x_ and x (which is x +
-        // some step)
-        system( x_, x1 );
-        system( x, x2 );
-        return ublas::norm_2( (x2 - x1)/dx_ );
-    }
-
-
-    /**
-     * @brief Makes the first guess. After this call the Newton method.
-     */
-    void correction_step(  )
-    {
-        // Get the jacobian at current point. Notice that in this method, we
-        // don't have to compute inverse of jacobian
-
-        vector_type direction( size_ );
-
-        // Now take the largest step possible such that the value of system at
-        // (x_ - step ) is lower than the value of system as x_.
-        vector_type nextState( size_ );
-
-        apply();
-
-        unsigned int i = 0;
-
-        double factor = 1e-2;
-        while( true )
-        {
-            i += 1;
-            compute_jacobians( x_, false );
-            // Make a move in either side of direction. In whichever direction
-            // the function decreases.
-            direction = ublas::prod( J_, f_ );
-            nextState = x_ - factor * direction;
-            if( ublas::norm_2( compute_at( nextState ) ) >= ublas::norm_2(compute_at(x_)))
-                factor = factor / 2.0;
-            else
-            {
-                cerr << "Correction term applied ";
-                x_ = nextState;
-                apply();
-                break;
-            }
-
-            if ( i > 20 )
-                break;
-        }
-    }
-
 public:
     const size_t size_;
 
diff --git a/moose-core/ksolve/SteadyStateBoost.cpp b/moose-core/ksolve/SteadyStateBoost.cpp
index 5d31165d66d52edc2d493aab558b429f151aff1b..98b671541edf292518cc1d70a90b38d936df6920 100644
--- a/moose-core/ksolve/SteadyStateBoost.cpp
+++ b/moose-core/ksolve/SteadyStateBoost.cpp
@@ -349,7 +349,7 @@ void SteadyState::setStoich( Id value )
     setupSSmatrix();
     double vol = LookupField< unsigned int, double >::get(stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
     pool_.setVolume( vol );
-    pool_.setStoich( stoichPtr, 0 );
+    pool_.setStoich( stoichPtr, nullptr );
     pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
     isInitialized_ = 1;
 }
@@ -561,6 +561,7 @@ void SteadyState::setupSSmatrix()
     total_.assign( nConsv, 0.0 );
 
     Id ksolve = Field< Id >::get( stoich_, "ksolve" );
+
     vector< double > nVec =
         LookupField< unsigned int, vector< double > >::get(
             ksolve,"nVec", 0 );
@@ -589,7 +590,7 @@ void SteadyState::classifyState( const double* T )
     double tot = 0.0;
     Stoich* s = reinterpret_cast< Stoich* >( stoich_.eref().data() );
     vector< double > nVec = LookupField< unsigned int, vector< double > >::get(
-                                s->getKsolve(), "nVec", 0 );
+            s->getKsolve(), "nVec", 0 );
     for ( unsigned int i = 0; i < numVarPools_; ++i )
         tot += nVec[i];
 
@@ -600,14 +601,14 @@ void SteadyState::classifyState( const double* T )
     for ( unsigned int i = 0; i < numVarPools_; ++i )
     {
         double orig = nVec[i];
-        if ( std::isnan( orig ) or std::isinf( orig ) )
+        if ( std::isnan( orig ) || std::isinf( orig ) )
         {
             cout << "Warning: SteadyState::classifyState: orig=nan\n";
             solutionStatus_ = 2; // Steady state OK, eig failed
             J.clear();
             return;
         }
-        if ( std::isnan( tot ) or std::isinf( tot ))
+        if ( std::isnan( tot ) || std::isinf( tot ))
         {
             cout << "Warning: SteadyState::classifyState: tot=nan\n";
             solutionStatus_ = 2; // Steady state OK, eig failed
@@ -622,73 +623,72 @@ void SteadyState::classifyState( const double* T )
         // Assign the rates for each mol.
         for ( unsigned int j = 0; j < numVarPools_; ++j )
         {
-            if( std::isnan( yprime[j] ) or std::isinf( yprime[j] ) )
+            if( std::isnan( yprime[j] ) || std::isinf( yprime[j] ) )
             {
-                cout << "Warning: Overflow/underflow. Can't continue " << endl;
                 solutionStatus_ = 2;
+                J.clear();
                 return;
             }
             J(i, j) = yprime[j];
         }
-    }
 
-    // Jacobian is now ready. Find eigenvalues.
-    ublas::vector< std::complex< double > > eigenVec ( J.size1() );
+        // Jacobian is now ready. Find eigenvalues.
+        ublas::vector< std::complex< double > > eigenVec ( J.size1() );
 
-    ublas::matrix< std::complex<double>, ublas::column_major >* vl, *vr;
-    vl = NULL; vr = NULL;
+        ublas::matrix< std::complex<double>, ublas::column_major >* vl, *vr;
+        vl = NULL; vr = NULL;
 
-    /*-----------------------------------------------------------------------------
-     *  INFO: Calling lapack routine geev to compute eigen vector of matrix J.
-     *
-     *  Argument 3 and 4 are left- and right-eigenvectors. Since we do not need
-     *  them, they are set to NULL. Argument 2 holds eigen-vector and result is
-     *  copied to it (output ).
-     *-----------------------------------------------------------------------------*/
-    int status = lapack::geev( J, eigenVec, vl, vr, lapack::optimal_workspace() );
+        /*-----------------------------------------------------------------------------
+         *  INFO: Calling lapack routine geev to compute eigen vector of matrix J.
+         *
+         *  Argument 3 and 4 are left- and right-eigenvectors. Since we do not need
+         *  them, they are set to NULL. Argument 2 holds eigen-vector and result is
+         *  copied to it (output ).
+         *-----------------------------------------------------------------------------*/
+        int status = lapack::geev( J, eigenVec, vl, vr, lapack::optimal_workspace() );
 
-    eigenvalues_.clear();
-    eigenvalues_.resize( numVarPools_, 0.0 );
-    if ( status != 0 )
-    {
-        cout << "Warning: SteadyState::classifyState failed to find eigenvalues. Status = " <<
-             status << endl;
-        solutionStatus_ = 2; // Steady state OK, eig classification failed
-    }
-    else     // Eigenvalues are ready. Classify state.
-    {
-        nNegEigenvalues_ = 0;
-        nPosEigenvalues_ = 0;
-        for ( unsigned int i = 0; i < numVarPools_; ++i )
+        eigenvalues_.clear();
+        eigenvalues_.resize( numVarPools_, 0.0 );
+        if ( status != 0 )
         {
-            std::complex<value_type> z = eigenVec[ i ];
-            double r = z.real();
-            nNegEigenvalues_ += ( r < -EPSILON );
-            nPosEigenvalues_ += ( r > EPSILON );
-            eigenvalues_[i] = r;
-            // We have a problem here because numVarPools_ usually > rank
-            // This means we have several zero eigenvalues.
+            cout << "Warning: SteadyState::classifyState failed to find eigenvalues. Status = " <<
+                status << endl;
+            solutionStatus_ = 2; // Steady state OK, eig classification failed
+        }
+        else     // Eigenvalues are ready. Classify state.
+        {
+            nNegEigenvalues_ = 0;
+            nPosEigenvalues_ = 0;
+            for ( unsigned int i = 0; i < numVarPools_; ++i )
+            {
+                std::complex<value_type> z = eigenVec[ i ];
+                double r = z.real();
+                nNegEigenvalues_ += ( r < -EPSILON );
+                nPosEigenvalues_ += ( r > EPSILON );
+                eigenvalues_[i] = r;
+                // We have a problem here because numVarPools_ usually > rank
+                // This means we have several zero eigenvalues.
+            }
+            if ( nNegEigenvalues_ == rank_ )
+                stateType_ = 0; // Stable
+            else if ( nPosEigenvalues_ == rank_ ) // Never see it.
+                stateType_ = 1; // Unstable
+            else  if (nPosEigenvalues_ == 1)
+                stateType_ = 2; // Saddle
+            else if ( nPosEigenvalues_ >= 2 )
+                stateType_ = 3; // putative oscillatory
+            else if ( nNegEigenvalues_ == ( rank_ - 1) && nPosEigenvalues_ == 0 )
+                stateType_ = 4; // one zero or unclassified eigenvalue. Messy.
+            else
+                stateType_ = 5; // Other
         }
-        if ( nNegEigenvalues_ == rank_ )
-            stateType_ = 0; // Stable
-        else if ( nPosEigenvalues_ == rank_ ) // Never see it.
-            stateType_ = 1; // Unstable
-        else  if (nPosEigenvalues_ == 1)
-            stateType_ = 2; // Saddle
-        else if ( nPosEigenvalues_ >= 2 )
-            stateType_ = 3; // putative oscillatory
-        else if ( nNegEigenvalues_ == ( rank_ - 1) && nPosEigenvalues_ == 0 )
-            stateType_ = 4; // one zero or unclassified eigenvalue. Messy.
-        else
-            stateType_ = 5; // Other
     }
 }
 
 static bool isSolutionValid( const vector< double >& x )
 {
-    for( size_t i = 0; i < x.size(); i++ )
+    for( auto &v : x )
     {
-        double v = x[i];
         if ( std::isnan( v ) || std::isinf( v ) )
         {
             cout << "Warning: SteadyState iteration gave nan/inf concs\n";
@@ -703,6 +703,19 @@ static bool isSolutionValid( const vector< double >& x )
     return true;
 }
 
+static bool isSolutionPositive( const vector< double >& x )
+{
+    for( auto &v : x )
+    {
+        if( v < 0.0 )
+        {
+            cout << "Warning: SteadyState iteration gave negative concs" << endl;
+            return false;
+        }
+    }
+    return true;
+}
+
 /**
  * The settle function computes the steady state nearest the initial
  * conditions.
@@ -773,7 +786,7 @@ void SteadyState::settle( bool forceSetup )
     int status = 1;
 
     // Find roots . If successful, set status to 0.
-    if( ss->find_roots_gnewton( ) )
+    if( ss->find_roots_gnewton( convergenceCriterion_, maxIter_ ) )
         status = 0;
 
     if ( status == 0 && isSolutionValid( ss->ri.nVec ) )
diff --git a/moose-core/ksolve/SteadyStateGsl.cpp b/moose-core/ksolve/SteadyStateGsl.cpp
index 5aa0941dacf666a66978c5b1ac66872dfe8731f9..eaec6d07ed115ed92fb4ec8cb9a4329a404412e6 100644
--- a/moose-core/ksolve/SteadyStateGsl.cpp
+++ b/moose-core/ksolve/SteadyStateGsl.cpp
@@ -383,7 +383,6 @@ void SteadyState::setStoich( Id value )
     pool_.setVolume( vol );
 
     pool_.setStoich( stoichPtr, nullptr );
-
     pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
     isInitialized_ = 1;
 }
@@ -936,13 +935,9 @@ int ss_func( const gsl_vector* x, void* params, gsl_vector* f )
     {
         double temp = op( gsl_vector_get( x, i ) );
         if ( isNaN( temp ) || isInfinity( temp ) )
-        {
             return GSL_ERANGE;
-        }
         else
-        {
             ri->nVec[i] = temp;
-        }
     }
     vector< double > vels;
     ri->pool->updateReacVelocities( &ri->nVec[0], vels );
diff --git a/moose-core/ksolve/VoxelPools.cpp b/moose-core/ksolve/VoxelPools.cpp
index 25e8ac7e875aed68b074c59d4b8a64dcb88a0593..53f5a1731621c21ded74935276779e2aa968c671 100644
--- a/moose-core/ksolve/VoxelPools.cpp
+++ b/moose-core/ksolve/VoxelPools.cpp
@@ -1,11 +1,11 @@
-/**********************************************************************
-** This program is part of 'MOOSE', the
-** Messaging Object Oriented Simulation Environment.
-**           Copyright (C) 2003-2010 Upinder S. Bhalla. and NCBS
-** It is made available under the terms of the
-** GNU Lesser General Public License version 2.1
-** See the file COPYING.LIB for the full notice.
-**********************************************************************/
+/*
+* This program is part of 'MOOSE', the
+* Messaging Object Oriented Simulation Environment.
+*           Copyright (C) 2003-2010 Upinder S. Bhalla. and NCBS
+* It is made available under the terms of the
+* GNU Lesser General Public License version 2.1
+* See the file COPYING.LIB for the full notice.
+*/
 
 #include "../basecode/header.h"
 
@@ -143,13 +143,20 @@ void VoxelPools::advance( const ProcInfo* p )
     // Variout stepper times are listed here:
     // https://www.boost.org/doc/libs/1_68_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.explicit_steppers
 
+    // Describe system to be used in boost solver calls.
     auto sys = [this](const vector_type_& dy, vector_type_& dydt, const double t) { 
         VoxelPools::evalRates(this, dy, dydt); };
 
-    // This is usually the default method for boost: Runge Kutta Fehlberg
-    if( method_ == "rk5")
-        odeint::integrate_const( rk_karp_stepper_type_()
-                , sys , Svec() , p->currTime - p->dt, p->currTime, p->dt);
+    // This is usually the default method. It works well in practice. Tested
+    // with steady-state solver. Closest to GSL rk5 .
+    if( method_ == "rk5" || method_ == "gsl" || method_ == "boost" )
+        odeint::integrate_const( 
+                make_dense_output( epsAbs_, epsRel_, odeint::runge_kutta_dopri5<vector_type_>() ) 
+                , sys , Svec() , p->currTime - p->dt , p->currTime , p->dt
+                );
+    else if( method_ == "rk5a" || method_ == "adaptive" )
+        odeint::integrate_adaptive( odeint::make_controlled<rk_dopri_stepper_type_>( epsAbs_, epsRel_ )
+                , sys , Svec() , p->currTime - p->dt , p->currTime, p->dt );
     else if( method_ == "rk2" )
         odeint::integrate_const( rk_midpoint_stepper_type_()
                 , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
@@ -203,17 +210,6 @@ int VoxelPools::gslFunc( double t, const double* y, double *dydt,
     VoxelPools* vp = reinterpret_cast< VoxelPools* >( params );
     // Stoich* s = reinterpret_cast< Stoich* >( params );
     double* q = const_cast< double* >( y ); // Assign the func portion.
-
-    // Assign the buffered pools
-    // Not possible because this is a static function
-    // Not needed because dydt = 0;
-    /*
-    double* b = q + s->getNumVarPools();
-    vector< double >::const_iterator sinit = Sinit_.begin() + s->getNumVarPools();
-    for ( unsigned int i = 0; i < s->getNumBufPools(); ++i )
-    	*b++ = *sinit++;
-    	*/
-
     vp->stoichPtr_->updateFuncs( q, t );
     vp->updateRates( y, dydt );
 #ifdef USE_GSL
diff --git a/moose-core/pymoose/CMakeLists.txt b/moose-core/pymoose/CMakeLists.txt
index d25a087569756a89e51281835928226820eebe5e..0a2a63885dbdac0e565d37df2e5fa49cdfa3f135 100644
--- a/moose-core/pymoose/CMakeLists.txt
+++ b/moose-core/pymoose/CMakeLists.txt
@@ -4,18 +4,18 @@ set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../cmake_modules/")
 
 find_package(PythonInterp REQUIRED)
 
-# Find Numpy 
+# Find Numpy
 find_package(NumPy REQUIRED)
 include_directories(${NUMPY_INCLUDE_DIRS})
 add_definitions( -std=c++11 )
 add_definitions(-DUSE_NUMPY)
 add_definitions(-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION)
 
-# set module extensiton. default is .so
+# set module extensiton. default is .so. Also check ../python/setup.cmake.py
 execute_process( COMMAND
     ${PYTHON_EXECUTABLE} -c
 "import importlib.machinery
-print(importlib.machinery.EXTENSION_SUFFIXES[0])"
+print(importlib.machinery.EXTENSION_SUFFIXES[-1])"
     OUTPUT_VARIABLE PYTHON_SO_EXTENSION
     OUTPUT_STRIP_TRAILING_WHITESPACE
     )
@@ -25,7 +25,7 @@ if(NOT PYTHON_SO_EXTENSION)
 endif()
 message(STATUS "Python so extension ${PYTHON_SO_EXTENSION}" )
 
-# TARGET 
+# TARGET
 set(PYMOOSE_SRCS
     moosemodule.cpp
     vec.cpp
@@ -63,7 +63,7 @@ else()
          COMPILE_DEFINITIONS "PYMOOSE"
          COMPILE_FLAGS "${PYTHON_INCLUDE_FLAGS}"
         )
-     
+
 endif()
 
 # Remove prefix lib from python module.
@@ -82,7 +82,7 @@ if(MACOSX AND ("${PYTHON_VERSION_MAJOR}" STREQUAL "2"))
 endif()
 
 # see issue #80
-if(HDF5_LIBRARY_DIRS)
+if(HDF5_FOUND AND WITH_NSDF)
     set_target_properties( _moose PROPERTIES LINK_FLAGS "-L${HDF5_LIBRARY_DIRS}" )
 endif()
 
diff --git a/moose-core/python/moose/moose.py b/moose-core/python/moose/moose.py
index b9acb706279e4744800bed8e1690419953f63c40..e6cf9bf7671de4bcf81dca56a7df495357dd09b5 100644
--- a/moose-core/python/moose/moose.py
+++ b/moose-core/python/moose/moose.py
@@ -29,10 +29,10 @@ try:
     import moose.neuroml2 as _neuroml2
 except Exception as e:
     nml2Import_ = False
-    nml2ImportError_ = '\n'.join( [ 
+    nml2ImportError_ = ' '.join( [ 
         "NML2 support is disabled because `libneuroml` and "
-        , "`pyneuroml` modules are not found."
-        , "     pip install pyneuroml libneuroml "
+        , "`pyneuroml` modules are not found.\n"
+        , "     $ pip install pyneuroml libneuroml \n"
         , " should fix it." 
         , " Actual error: %s " % e ]
         )
@@ -229,27 +229,23 @@ def mergeChemModel(src, des):
 
 # NML2 reader and writer function.
 def mooseReadNML2( modelpath ):
-    """Read NeuroML model (version 2).
-
+    """Read NeuroML model (version 2) and return reader object.
     """
     global nml2Import_
-    if nml2Import_:
-        reader = _neuroml2.NML2Reader( )
-        reader.read( modelpath )
-        return reader
-    else:
-        mu.info( nml2ImportError_ )
-        mu.warn( "Could not load NML2 support. Doing nothing" )
-        return False
+    if not nml2Import_:
+        mu.warn( nml2ImportError_ )
+        raise RuntimeError( "Could not load NML2 support." )
+
+    reader = _neuroml2.NML2Reader( )
+    reader.read( modelpath )
+    return reader
 
 def mooseWriteNML2( outfile ):
-    mu.warn( "Writing to NML2 is not supported yet" )
+    raise NotImplementedError( "Writing to NML2 is not supported yet" )
 
 ################################################################
 # Wrappers for global functions
 ################################################################
-
-
 def pwe():
     """Print present working element. Convenience function for GENESIS
     users. If you want to retrieve the element in stead of printing
diff --git a/moose-core/python/moose/neuroml/ChannelML.py b/moose-core/python/moose/neuroml/ChannelML.py
index 60a566fdfb7fa6edd5b51dbb6d10e20081d6117a..cc2c44437609d793bb0634f75632db66b0127c54 100644
--- a/moose-core/python/moose/neuroml/ChannelML.py
+++ b/moose-core/python/moose/neuroml/ChannelML.py
@@ -15,14 +15,13 @@ readChannelML(...) / readSynapseML to load from an xml.etree xml element (could
 
 from __future__ import print_function
 from xml.etree import ElementTree as ET
-import string
-import os, sys
+import os
+import sys
 import math
-
 import moose
 from moose.neuroml import utils
-from moose import utils as moose_utils
-from moose import print_utils as pu
+import moose.utils as mu
+import moose.print_utils as pu
 
 class ChannelML():
 
@@ -497,7 +496,7 @@ def make_new_synapse(syn_name, postcomp, syn_name_full, nml_params):
     ## connect the SimpleSynHandler or the STDPSynHandler to the SynChan (double exp)
     moose.connect( synhandler, 'activationOut', syn, 'activation' )
     # mgblock connections if required
-    childmgblock = moose_utils.get_child_Mstring(syn,'mgblockStr')
+    childmgblock = mu.get_child_Mstring(syn,'mgblockStr')
     #### connect the post compartment to the synapse
     if childmgblock.value=='True': # If NMDA synapse based on mgblock, connect to mgblock
         mgblock = moose.Mg_block(syn.path+'/mgblock')
diff --git a/moose-core/python/moose/neuroml/MorphML.py b/moose-core/python/moose/neuroml/MorphML.py
index 0a5d0c156a25fc7977a6fa584f15797589130431..d34a0c12c305276fe7e4c7b4d3d6907754773293 100644
--- a/moose-core/python/moose/neuroml/MorphML.py
+++ b/moose-core/python/moose/neuroml/MorphML.py
@@ -259,15 +259,18 @@ class MorphML():
                 if cableid in self.intFireCableIds:
                     mechanismname = self.intFireCableIds[cableid]
                 if mechanismname is not None: # this cableid is an intfire
-                    ## create LIF (subclass of Compartment) and set to default values
+                    # create LIF (subclass of Compartment) and set to default values
                     moosecomp = moose.LIF(moosecomppath)
-                    moosechannel = moose.Neutral('/library/'+mechanismname)
+                    mname = '/library/' + mechanismname
+                    moosechannel = moose.element(mname) if moose.exists(mname) else moose.Neutral(mname)
+                    # Mstring values are 'string'; make sure to convert them to
+                    # float else it will seg-fault with python3+
                     moosechannelval = moose.Mstring(moosechannel.path+'/vReset')
-                    moosecomp.vReset = moosechannelval.value
+                    moosecomp.vReset = float(moosechannelval.value)
                     moosechannelval = moose.Mstring(moosechannel.path+'/thresh')
-                    moosecomp.thresh = moosechannelval.value
+                    moosecomp.thresh = float( moosechannelval.value )
                     moosechannelval = moose.Mstring(moosechannel.path+'/refracT')
-                    moosecomp.refractoryPeriod = moosechannelval.value
+                    moosecomp.refractoryPeriod = eval(moosechannelval.value)
                     ## refracG is currently not supported by moose.LIF
                     ## when you implement it, check if refracG or g_refrac
                     ## is a conductance density or a conductance, I think the former
diff --git a/moose-core/python/moose/neuroml/NetworkML.py b/moose-core/python/moose/neuroml/NetworkML.py
index d7f30692eb11115bf00e9d2fdd59d760b16c7e8b..14ae5f53598e34a5ce1ef96ea44d74f83bcaf2eb 100644
--- a/moose-core/python/moose/neuroml/NetworkML.py
+++ b/moose-core/python/moose/neuroml/NetworkML.py
@@ -179,7 +179,7 @@ class NetworkML():
             model_filenames = (cellname+'.xml', cellname+'.morph.xml')
             success = False
             for model_filename in model_filenames:
-                model_path = find_first_file(model_filename,self.model_dir)
+                model_path = find_first_file(model_filename, self.model_dir)
                 if model_path is not None:
                     cellDict = mmlR.readMorphMLFromFile(model_path, self.params)
                     success = True
diff --git a/moose-core/python/moose/neuroml/NeuroML.py b/moose-core/python/moose/neuroml/NeuroML.py
index ceddf21d811b33e77f2828a269a3942b15c65085..5d9d8edf4799f6947ca1e97d6720d525e59528b3 100644
--- a/moose-core/python/moose/neuroml/NeuroML.py
+++ b/moose-core/python/moose/neuroml/NeuroML.py
@@ -1,8 +1,12 @@
 # -*- coding: utf-8 -*-
-## Description: class NeuroML for loading NeuroML from single file into MOOSE
-## Version 1.0 by Aditya Gilra, NCBS, Bangalore, India, 2011 for serial MOOSE
-## Version 1.5 by Niraj Dudani, NCBS, Bangalore, India, 2012, ported to parallel MOOSE
-## Version 1.6 by Aditya Gilra, NCBS, Bangalore, India, 2012, further changes for parallel MOOSE
+from __future__ import print_function, division
+
+# Description: class NeuroML for loading NeuroML from single file into MOOSE
+# Version 1.0 by Aditya Gilra, NCBS, Bangalore, India, 2011 for serial MOOSE
+# Version 1.5 by Niraj Dudani, NCBS, Bangalore, India, 2012, ported to parallel MOOSE
+# Version 1.6 by Aditya Gilra, NCBS, Bangalore, India, 2012, further changes for parallel MOOSE
+# Maintainer: Dilawar Singh, dilawars@ncbs.res.in 
+#  - python3 related fixes.
 
 """
 NeuroML.py is the preferred interface to read NeuroML files.
@@ -15,7 +19,8 @@ readNeuroMLFromFile(...) to load NeuroML from a file:
 
 OR
 
-(b) the file could have only L3 (network) with L2 (channels/synapses) and L1 (cells) spread over multiple files;
+(b) the file could have only L3 (network) with L2 (channels/synapses) and L1
+(cells) spread over multiple files;
  these multiple files should be in the same or lower-level directory
  named as <chan/syn_name>.xml or <cell_name>.xml or <cell_name>.morph.xml
  (essentially as generated by neuroConstruct's export).
@@ -44,7 +49,6 @@ In [2]: import moose.neuroml
 In [3]: moose.neuroml.loadNeuroML_L123('Generated.net.xml')
 """
 
-from __future__ import print_function
 import moose
 from moose.utils import *
 from xml.etree import ElementTree as ET
@@ -57,32 +61,27 @@ from moose.neuroml.utils import *
 import sys
 from os import path
 
-import logging
-console = logging.StreamHandler()
-console.setLevel(logging.INFO)
-formatter = logging.Formatter('%(name)-12s: %(levelname)-8s %(message)s')
-console.setFormatter(formatter)
-_logger = logging.getLogger('moose.nml.neuroml')
-_logger.addHandler(console)
-
-
 class NeuroML():
 
     def __init__(self):
         pass
 
-    def readNeuroMLFromFile(self,filename,params={},cellsDict={}):
+    def readNeuroMLFromFile(self, filename, params={}, cellsDict={}):
         """
         For the format of params required to tweak what cells are loaded,
          refer to the doc string of NetworkML.readNetworkMLFromFile().
         Returns (populationDict,projectionDict),
          see doc string of NetworkML.readNetworkML() for details.
         """
-        _logger.info("Loading neuroml file %s " % filename)
+        mu.info("Loading neuroml file %s " % filename)
         moose.Neutral('/library') # creates /library in MOOSE tree; elif present, wraps
         tree = ET.parse(filename)
         root_element = tree.getroot()
-        self.model_dir = path.dirname( path.abspath( filename ) )
+
+        # if model_path is given in params, use it else use the directory of NML
+        # as model_dir.
+        self.model_dir = params.get('model_dir', path.dirname(path.abspath(filename)))
+
         if 'lengthUnits' in list(root_element.attrib.keys()):
             self.lengthUnits = root_element.attrib['lengthUnits']
         else:
@@ -110,14 +109,14 @@ class NeuroML():
                         self.temperature = float(tag.find('.//{'+meta_ns+'}value').text)
                         self.temperature_default = False
         if self.temperature_default:
-            _logger.info("Using default temperature of %s degree Celsius" % self.temperature)
+            mu.info("Using default temperature of %s degree Celsius" % self.temperature)
         self.nml_params = {
                 'temperature':self.temperature,
                 'model_dir':self.model_dir,
         }
 
 
-        _logger.debug("Loading channels and synapses into MOOSE /library ...")
+        mu.debug("Loading channels and synapses into MOOSE /library ...")
         cmlR = ChannelML(self.nml_params)
         for channels in root_element.findall('.//{'+neuroml_ns+'}channels'):
             self.channelUnits = channels.attrib['units']
@@ -131,7 +130,7 @@ class NeuroML():
             for ionConc in channels.findall('.//{'+cml_ns+'}ion_concentration'):
                 cmlR.readIonConcML(ionConc,units=self.channelUnits)
 
-        _logger.debug("Loading cell definitions into MOOSE /library ...")
+        mu.debug("Loading cell definitions into MOOSE /library ...")
         mmlR = MorphML(self.nml_params)
         self.cellsDict = cellsDict
         for cells in root_element.findall('.//{'+neuroml_ns+'}cells'):
@@ -145,7 +144,7 @@ class NeuroML():
             and root_element.find('.//{'+nml_ns+'}populations') is None:
             return (self.cellsDict,'no populations (L3 NetworkML) found.')
         else:
-            _logger.debug("Loading individual cells into MOOSE root ... ")
+            mu.debug("Loading individual cells into MOOSE root ... ")
             nmlR = NetworkML(self.nml_params)
             return nmlR.readNetworkML(root_element,self.cellsDict,\
                     params=params,lengthUnits=self.lengthUnits)
@@ -163,6 +162,6 @@ def loadNeuroML_L123(filename):
 
 if __name__ == "__main__":
     if len(sys.argv)<2:
-        _logger.error("You need to specify the neuroml filename.")
+        mu.error("You need to specify the neuroml filename.")
         sys.exit(1)
     print(loadNeuroML_L123(sys.argv[1]))
diff --git a/moose-core/python/moose/neuroml2/reader.py b/moose-core/python/moose/neuroml2/reader.py
index 1fdf4fc197ae8976016156cc632a8bee5c53b639..75c9a7f277d67889bb381180f4fa4044e4fbc0a7 100644
--- a/moose-core/python/moose/neuroml2/reader.py
+++ b/moose-core/python/moose/neuroml2/reader.py
@@ -9,48 +9,11 @@
 # Version: 
 # Last-Updated: 15 Jan 2018, pgleeson
 #               16 Jan 2018, dilawar, python3 compatible imports.
-#
-# URL:
-# Keywords:
-# Compatibility:
-#
-
-# Commentary:
-#
-#
-#
-#
-
-# Change log:
-#
-#
-#
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License as
-# published by the Free Software Foundation; either version 3, or
-# (at your option) any later version.
-# 
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-# General Public License for more details.
-# 
-# You should have received a copy of the GNU General Public License
-# along with this program; see the file COPYING.  If not, write to
-# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
-# Floor, Boston, MA 02110-1301, USA.
-#
-#
-
-# Code:
-"""Implementation of reader for NeuroML 2 models.
-
+#             
 
+"""Implementation of reader for NeuroML 2 models.
 TODO: handle morphologies of more than one segment...
-
 """
-
 from __future__ import print_function, division, absolute_import
 
 try:
@@ -64,13 +27,19 @@ import math
 import numpy as np
 import neuroml as nml
 from pyneuroml import pynml
-
 import moose
 import moose.utils as mu
 
 from .units import SI
 from . import hhfit
 
+def _unique( ls ):
+    res = [ ]
+    for l in ls:
+        if l not in res:
+            res.append( l )
+    return res
+
 def sarea(comp):
     """
     Return the surface area of compartment from length and
@@ -116,15 +85,13 @@ def setEk(comp, erev):
 def getSegments(nmlcell, component, sg_to_segments):
     """Get the list of segments the `component` is applied to"""
     sg = component.segment_groups
-    #seg = component.segment
     if sg is None:
-            segments = nmlcell.morphology.segments
-    elif sg == 'all': # Special case
+        segments = nmlcell.morphology.segments
+    elif sg == 'all': 
         segments = [seg for seglist in sg_to_segments.values() for seg in seglist]
     else:
         segments = sg_to_segments[sg]
-        
-    return list(set(segments))
+    return _unique(segments)
 
 
 class NML2Reader(object):
@@ -148,10 +115,10 @@ class NML2Reader(object):
         self.filename = None        
         self.nml_to_moose = {} # NeuroML object to MOOSE object
         self.moose_to_nml = {} # Moose object to NeuroML object
-        self.proto_cells = {} # map id to prototype cell in moose
-        self.proto_chans = {} # map id to prototype channels in moose
-        self.proto_pools = {} # map id to prototype pools (Ca2+, Mg2+)
-        self.includes = {} # Included files mapped to other readers
+        self.proto_cells = {}  # map id to prototype cell in moose
+        self.proto_chans = {}  # map id to prototype channels in moose
+        self.proto_pools = {}  # map id to prototype pools (Ca2+, Mg2+)
+        self.includes = {}     # Included files mapped to other readers
         self.lib = moose.Neutral('/library')
         self.id_to_ionChannel = {}
         self._cell_to_sg = {} # nml cell to dict - the dict maps segment groups to segments
@@ -198,11 +165,14 @@ class NML2Reader(object):
         return self.cells_in_populations[pop_id][index]
     
     def getComp(self, pop_id, cellIndex, segId):
-        return moose.element('%s/%s/%s/%s' % (self.lib.path, pop_id, cellIndex, self.seg_id_to_comp_name[self.pop_to_cell_type[pop_id]][segId]))
+        return moose.element('%s/%s/%s/%s' % ( self.lib.path, pop_id, cellIndex
+            , self.seg_id_to_comp_name[self.pop_to_cell_type[pop_id]][segId])
+            )
             
     def createPopulations(self):
         for pop in self.network.populations:
-            mpop = moose.Neutral('%s/%s' % (self.lib.path, pop.id))
+            ep = '%s/%s' % (self.lib.path, pop.id)
+            mpop = moose.element(ep) if moose.exists(ep) else moose.Neutral(ep)
             self.cells_in_populations[pop.id] ={}
             for i in range(pop.size):
                 mu.info("Creating %s/%s instances of %s under %s"%(i,pop.size,pop.component, mpop))
@@ -228,12 +198,15 @@ class NML2Reader(object):
         for il in self.network.input_lists:
             for ii in il.input:
                 input = self.getInput(il.component)
-                moose.connect(input, 'output', self.getComp(il.populations,ii.get_target_cell_id(),ii.get_segment_id()), 'injectMsg')
+                moose.connect(input, 'output'
+                        , self.getComp(il.populations,ii.get_target_cell__hash(),ii.get_segment__hash())
+                        , 'injectMsg')
             
 
     def createCellPrototype(self, cell, symmetric=True):
         """To be completed - create the morphology, channels in prototype"""
-        nrn = moose.Neuron('%s/%s' % (self.lib.path, cell.id))
+        ep = '%s/%s' % (self.lib.path, cell.id)
+        nrn = moose.element(ep) if moose.exists(ep) else moose.Neuron(ep)
         self.proto_cells[cell.id] = nrn
         self.nml_to_moose[cell] = nrn
         self.moose_to_nml[nrn] = cell
@@ -262,11 +235,13 @@ class NML2Reader(object):
         self.seg_id_to_comp_name[nmlcell.id]={}
         for seg in segments:
             if seg.name is not None:
-                id_to_comp[seg.id] = compclass('%s/%s' % (cellpath, seg.name))
+                ep = '%s/%s' % (cellpath, seg.name)
+                id_to_comp[seg.id] = moose.element(ep) if moose.exists(ep) else compclass(ep)
                 self.seg_id_to_comp_name[nmlcell.id][seg.id] = seg.name
             else:
                 name = 'comp_%s'%seg.id
-                id_to_comp[seg.id] = compclass('%s/%s' % (cellpath, name))
+                ep = '%s/%s' % (cellpath, name)
+                id_to_comp[seg.id] = moose.element(ep) if moose.exists(ep) else compclass(ep)
                 self.seg_id_to_comp_name[nmlcell.id][seg.id] = name
         # Now assign the positions and connect up axial resistance
         if not symmetric:
@@ -286,13 +261,15 @@ class NML2Reader(object):
                 if parent:
                     p0 = parent.distal
                 else:
-                    raise Exception('No proximal point and no parent segment for segment: name=%s, id=%s' % (segment.name, segment.id))
-            comp.x0, comp.y0, comp.z0 = (x * self.lunit for x in map(float, (p0.x, p0.y, p0.z)))
+                    raise Exception(
+                    'No proximal point and no parent segment for segment:'+\
+                            'name=%s, id=%s' % (segment.name, segment.id)
+                        )
+            comp.x0, comp.y0, comp.z0 = (x*self.lunit for x in map(float, (p0.x, p0.y, p0.z)))
             p1 = segment.distal
-            comp.x, comp.y, comp.z = (x * self.lunit for x in map(float, (p1.x, p1.y, p1.z)))
-            comp.length = np.sqrt((comp.x - comp.x0)**2
-                                  + (comp.y - comp.y0)**2
-                                  + (comp.z - comp.z0)**2)
+            comp.x, comp.y, comp.z = (x*self.lunit for x in map(float, (p1.x, p1.y, p1.z)))
+            comp.length = np.sqrt((comp.x-comp.x0)**2+(comp.y-comp.y0)**2+(comp.z-comp.z0)**2)
+
             # This can pose problem with moose where both ends of
             # compartment have same diameter. We are averaging the two
             # - may be splitting the compartment into two is better?
@@ -377,10 +354,15 @@ class NML2Reader(object):
                     proto_pool = innerReader.proto_pools[species.concentrationModel]
                     break
         if not proto_pool:
-            raise Exception('No prototype pool for %s referred to by %s' % (species.concentration_model, species.id))
+            raise Exception('No prototype pool for %s referred to by %s' % ( 
+                    species.concentration_model, species.id)
+                )
         pool_id = moose.copy(proto_pool, comp, species.id)
         pool = moose.element(pool_id)
-        pool.B = pool.B / (np.pi * compartment.length * (0.5 * compartment.diameter + pool.thickness) * (0.5 * compartment.diameter - pool.thickness))        
+        pool.B = pool.B / (np.pi * compartment.length * ( 
+            0.5 * compartment.diameter + pool.thickness) * 
+            (0.5 * compartment.diameter - pool.thickness)
+            )
         return pool
 
     def importAxialResistance(self, nmlcell, intracellularProperties):
@@ -417,8 +399,11 @@ class NML2Reader(object):
                     mu.info("Using %s to evaluate rate"%ct.name)
                     rate = []
                     for v in tab:
-                        vals = pynml.evaluate_component(ct,req_variables={'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()})
-                        '''mu.info vals'''
+                        vals = pynml.evaluate_component(ct
+                                , req_variables = 
+                                    {'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()}
+                                )
+                        # mu.info vals
                         if 'x' in vals:
                             rate.append(vals['x'])
                         if 't' in vals:
@@ -483,32 +468,6 @@ class NML2Reader(object):
         moose.connect(chan, 'channel', comp, 'channel')
         return chan    
 
-    '''
-    def importIncludes(self, doc):        
-        for include in doc.include:
-            if self.verbose:
-                mu.info(self.filename, 'Loading include', include)
-            error = None
-            inner = NML2Reader(self.verbose)
-            paths = [include.href, os.path.join(os.path.dirname(self.filename), include.href)]
-            for path in paths:
-                try:
-                    inner.read(path)                    
-                    if self.verbose:
-                        mu.info(self.filename, 'Loaded', path, '... OK')
-                except IOError as e:
-                    error = e
-                else:
-                    self.includes[include.href] = inner
-                    self.id_to_ionChannel.update(inner.id_to_ionChannel)
-                    self.nml_to_moose.update(inner.nml_to_moose)
-                    self.moose_to_nml.update(inner.moose_to_nml)
-                    error = None
-                    break
-            if error:
-                mu.info(self.filename, 'Last exception:', error)
-                raise IOError('Could not read any of the locations: %s' % (paths))'''
-                
     def _is_standard_nml_rate(self, rate):
         return rate.type=='HHExpLinearRate' \
                or rate.type=='HHExpRate' or \
@@ -517,7 +476,7 @@ class NML2Reader(object):
 
     def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000):
         mchan = moose.HHChannel('%s/%s' % (self.lib.path, chan.id))
-        mgates = map(moose.element, (mchan.gateX, mchan.gateY, mchan.gateZ))
+        mgates = [moose.element(x) for x in [mchan.gateX, mchan.gateY, mchan.gateZ]]
         assert(len(chan.gate_hh_rates) <= 3) # We handle only up to 3 gates in HHCHannel
         
         if self.verbose:
@@ -551,19 +510,25 @@ class NML2Reader(object):
                 if ngate.q10_settings.type == 'q10Fixed':
                     q10_scale= float(ngate.q10_settings.fixed_q10)
                 elif ngate.q10_settings.type == 'q10ExpTemp':
-                    q10_scale = math.pow(float(ngate.q10_settings.q10_factor),(self._getTemperature()- SI(ngate.q10_settings.experimental_temp))/10)
-                    #mu.info('Q10: %s; %s; %s; %s'%(ngate.q10_settings.q10_factor, self._getTemperature(),SI(ngate.q10_settings.experimental_temp),q10_scale))
+                    q10_scale = math.pow( float(ngate.q10_settings.q10_factor)
+                            , (self._getTemperature()- SI(ngate.q10_settings.experimental_temp))/10
+                            )
                 else:
-                    raise Exception('Unknown Q10 scaling type %s: %s'%(ngate.q10_settings.type,ngate.q10_settings))
+                    raise Exception('Unknown Q10 scaling type %s: %s'%( 
+                        ngate.q10_settings.type,ngate.q10_settings)
+                        )
                     
             if self.verbose:
-                mu.info(' === Gate: %s; %s; %s; %s; %s; scale=%s'%(ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale))
+                mu.info(' === Gate: %s; %s; %s; %s; %s; scale=%s'% ( 
+                    ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale)
+                    )
                 
             if (fwd is not None) and (rev is not None):
                 alpha = self.calculateRateFn(fwd, vmin, vmax, vdivs)
                 beta = self.calculateRateFn(rev, vmin, vmax, vdivs)
                 mgate.tableA = q10_scale * (alpha)
                 mgate.tableB = q10_scale * (alpha + beta)
+
             # Assuming the meaning of the elements in GateHHTauInf ...
             if hasattr(ngate,'time_course') and hasattr(ngate,'steady_state') \
                and (ngate.time_course is not None) and (ngate.steady_state is not None):
@@ -587,16 +552,25 @@ class NML2Reader(object):
         return mchan
 
     def createPassiveChannel(self, chan):
-        mchan = moose.Leakage('%s/%s' % (self.lib.path, chan.id))
+        epath = '%s/%s' % (self.lib.path, chan.id)
+        if moose.exists( epath ):
+            mchan = moose.element(epath)
+        else:
+            mchan = moose.Leakage(epath)
         if self.verbose:
             mu.info(self.filename, 'Created', mchan.path, 'for', chan.id)
         return mchan
 
     def importInputs(self, doc):
-        minputs = moose.Neutral('%s/inputs' % (self.lib.path))
-        for pg_nml in doc.pulse_generators:
+        epath = '%s/inputs' % (self.lib.path)
+        if moose.exists( epath ):
+            minputs = moose.element( epath )
+        else:
+            minputs = moose.Neutral( epath )
 
-            pg = moose.PulseGen('%s/%s' % (minputs.path, pg_nml.id))
+        for pg_nml in doc.pulse_generators:
+            epath = '%s/%s' % (minputs.path, pg_nml.id)
+            pg = moose.element(epath) if moose.exists(epath) else moose.PulseGen(epath)
             pg.firstDelay = SI(pg_nml.delay)
             pg.firstWidth = SI(pg_nml.duration)
             pg.firstLevel = SI(pg_nml.amplitude)
@@ -639,14 +613,10 @@ class NML2Reader(object):
         ca.CaBasal = SI(concModel.restingConc)
         ca.tau = SI(concModel.decayConstant)
         ca.thick = SI(concModel.shellThickness)
-        ca.B = 5.2e-6 # B = 5.2e-6/(Ad) where A is the area of the shell and d is thickness - must divide by shell volume when copying
+        ca.B = 5.2e-6 # B = 5.2e-6/(Ad) where A is the area of the 
+                      # shell and d is thickness - must divide by 
+                      # shell volume when copying
         self.proto_pools[concModel.id] = ca
         self.nml_to_moose[concModel.id] = ca
         self.moose_to_nml[ca] = concModel
         logger.debug('Created moose element: %s for nml conc %s' % (ca.path, concModel.id))
-
-
-
-
-# 
-# reader.py ends here
diff --git a/moose-core/python/setup.cmake.py b/moose-core/python/setup.cmake.py
index e6bbd1a7f749780f16c11b508e8ec0dc4dfbe183..c7254b402121afb4ebe16364f81c2c3204914599 100644
--- a/moose-core/python/setup.cmake.py
+++ b/moose-core/python/setup.cmake.py
@@ -32,11 +32,12 @@ with open( os.path.join( script_dir, 'VERSION'), 'r' ) as f:
     version = f.read( )
 print( 'Got %s from VERSION file' % version )
 
-# importlib is available only for python3.
+# importlib is available only for python3. Since we build wheels, prefer .so
+# extension. This way a wheel built by any python3.x will work with any python3.
 suffix = '.so'
 try:
     import importlib.machinery
-    suffix = importlib.machinery.EXTENSION_SUFFIXES[0]
+    suffix = importlib.machinery.EXTENSION_SUFFIXES[-1]
 except Exception as e:
     print( '[WARN] Failed to determine importlib suffix' )
     suffix = '.so'
@@ -54,9 +55,9 @@ setup(
         packages               = [ 'rdesigneur', 'moose'
                                     , 'moose.SBML', 'moose.genesis'
                                     , 'moose.neuroml', 'moose.neuroml2'
-                                    , 'moose.chemUtil', 'moose.chemMerge' 
+                                    , 'moose.chemUtil', 'moose.chemMerge'
                                 ],
-        install_requires       = [ 'python-libsbml', 'numpy' ],
+        install_requires       = [ 'numpy' ],
         package_dir            = { 'moose' : 'moose', 'rdesigneur' : 'rdesigneur' },
         package_data           = { 'moose' : ['_moose' + suffix, 'neuroml2/schema/NeuroMLCoreDimensions.xml'] },
         )
diff --git a/moose-core/tests/python/test_Xreacs2.py b/moose-core/tests/python/test_Xreacs2.py
index 4089501cc4b8fcfbc27a8484ae5b2d9e7f6be87a..2fb36cfb16f05a3e4669b6f201f30867330a25f1 100644
--- a/moose-core/tests/python/test_Xreacs2.py
+++ b/moose-core/tests/python/test_Xreacs2.py
@@ -7,7 +7,8 @@ import moose.fixXreacs as fixXreacs
 
 def countCrossings( plot, thresh ):
     vec = moose.element( plot ).vector
-    #print (vec[:-1] < thresh)
+    print( vec )
+    #  print (vec[:-1] <= thresh)
     return sum( (vec[:-1] < thresh) * (vec[1:] >= thresh ) )
 
 def main( standalone = False ):
@@ -43,7 +44,8 @@ def main( standalone = False ):
     moose.reinit()
     moose.start( runtime )
     # I don't have an analytic way to assess oscillations
-    assert( countCrossings( '/model/graphs/conc2/M.Co', 0.001 ) == 4 )
+    nCrossings = countCrossings( '/model/graphs/conc2/M.Co', 0.001 )
+    assert( nCrossings == 4 ), "Expected 4, got %d" % nCrossings
     moose.delete( '/model' )
 
 # Run the 'main' if this script is executed standalone.
diff --git a/moose-core/tests/python/test_dose_response.py b/moose-core/tests/python/test_dose_response.py
new file mode 100644
index 0000000000000000000000000000000000000000..4a8973013f1c1ed7b3c220bfb10f626c642fb8e7
--- /dev/null
+++ b/moose-core/tests/python/test_dose_response.py
@@ -0,0 +1,116 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function, division
+
+## Makes and plots the dose response curve for bistable models
+## Author: Sahil Moza
+## June 26, 2014
+## Update:
+## Friday 14 September 2018 05:48:42 PM IST
+## Tunrned into a light-weight test by Dilawar Singh 
+
+import os
+import sys
+import moose
+import numpy as np
+
+sdir_ = os.path.dirname( os.path.realpath( __file__ ) )
+vals_ = [ ]
+
+def setupSteadyState(simdt,plotDt):
+
+    ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
+    stoich = moose.Stoich( '/model/kinetics/stoich' )
+    stoich.compartment = moose.element('/model/kinetics')
+    stoich.ksolve = ksolve
+    stoich.path = "/model/kinetics/##"
+    state = moose.SteadyState( '/model/kinetics/state' )
+    moose.reinit()
+    state.stoich = stoich
+    state.showMatrices()
+    state.convergenceCriterion = 1e-8
+    return ksolve, state
+
+def parseModelName(fileName):
+    pos1=fileName.rfind('/')
+    pos2=fileName.rfind('.')
+    directory=fileName[:pos1]
+    prefix=fileName[pos1+1:pos2]
+    suffix=fileName[pos2+1:len(fileName)]
+    return directory, prefix, suffix
+
+# Solve for the steady state
+def getState( ksolve, state, vol):
+    scale = 1.0 / ( vol * 6.022e23 )
+    moose.reinit()
+    state.randomInit() # Removing random initial condition to systematically make Dose reponse curves.
+    moose.start( 2.0 ) # Run the model for 2 seconds.
+    a = moose.element( '/model/kinetics/a' ).conc
+    vals_.append( a )
+    state.settle()
+    
+    vector = []
+    a = moose.element( '/model/kinetics/a' ).conc
+    for x in ksolve.nVec[0]:
+        vector.append( x * scale)
+    failedSteadyState = any([np.isnan(x) for x in vector])
+    if not (failedSteadyState):
+        return state.stateType, state.solutionStatus, a, vector
+
+
+def main():
+    # Setup parameters for simulation and plotting
+    simdt= 1e-2
+    plotDt= 1
+
+    # Factors to change in the dose concentration in log scale
+    factorExponent = 10  ## Base: ten raised to some power.
+    factorBegin = -10
+    factorEnd = 11
+    factorStepsize = 2
+    factorScale = 10.0 ## To scale up or down the factors
+
+    # Load Model and set up the steady state solver.
+    # model = sys.argv[1] # To load model from a file.
+    model = os.path.join( sdir_, 'chem_models/19085.cspace' )
+    modelPath, modelName, modelType = parseModelName(model)
+    outputDir = modelPath
+    
+    modelId = moose.loadModel(model, 'model', 'ee')
+    dosePath = '/model/kinetics/b/DabX' # The dose entity
+
+    ksolve, state = setupSteadyState( simdt, plotDt)
+    vol = moose.element( '/model/kinetics' ).volume
+    iterInit = 100
+    solutionVector = []
+    factorArr = []
+    
+    enz = moose.element(dosePath)
+    init = float(enz.kcat) # Dose parameter
+    
+    # Change Dose here to .
+    for factor in range(factorBegin, factorEnd, factorStepsize ):
+        scale = factorExponent ** (factor/factorScale) 
+        enz.kcat = init * scale     
+        print( "scale={:.3f}\tkcat={:.3f}".format( scale, enz.kcat) )
+        for num in range(iterInit):
+            stateType, solStatus, a, vector = getState( ksolve, state, vol)
+            if solStatus == 0:
+                #solutionVector.append(vector[0]/sum(vector))
+                solutionVector.append(a)
+                factorArr.append(scale)   
+
+    expected = (0.001411, 0.00092559)
+    got = ( np.mean(vals_), np.std(vals_) )
+    assert np.isclose(got, expected, atol=1e-4).all(), "Got %s, expected %s" % (got, expected)
+    print( '[TEST1] Passed. Concentration of a is same' )
+                
+    joint = np.array([factorArr, solutionVector])
+    joint = joint[:,joint[1,:].argsort()]
+    got = np.mean( joint ), np.std( joint )
+    expected = (1.2247, 2.46)
+    # Close upto 2 decimal place is good enough.
+    assert np.isclose(got, expected, atol=1e-2).all(), "Got %s, expected %s" % (got, expected)
+    print( joint )
+
+if __name__ == '__main__':
+    main()
diff --git a/moose-core/tests/python/test_neuroml2.py b/moose-core/tests/python/test_neuroml2.py
index 8460c092195f22de18a256cc7371045a8e0a192c..541d75806ae342d0129c9fe1637223bbf1776d73 100644
--- a/moose-core/tests/python/test_neuroml2.py
+++ b/moose-core/tests/python/test_neuroml2.py
@@ -14,6 +14,18 @@ import numpy as np
 
 SCRIPT_DIR = os.path.dirname( os.path.realpath( __file__ ) )
 
+# check if neuroml working properly.
+# NOTE: This script does not work with python3 
+# See https://github.com/NeuroML/NeuroML2/issues/116 . If this bug is fixed then
+# remove this code block.
+import neuroml as nml
+a = nml.nml.nml.IonChannel()
+try:
+    b = {a : 1 }
+except TypeError as e:
+    print( 'Failed due to https://github.com/NeuroML/NeuroML2/issues/116' ) 
+    quit( 0 )
+
 def run( nogui = True ):
     global SCRIPT_DIR
     filename = os.path.join(SCRIPT_DIR, 'test_files/passiveCell.nml' )
diff --git a/moose-core/tests/python/test_steady_state_solver.py b/moose-core/tests/python/test_steady_state_solver.py
new file mode 100644
index 0000000000000000000000000000000000000000..fd3e13e1b611e9e84bde2d53e3de56585785a9ea
--- /dev/null
+++ b/moose-core/tests/python/test_steady_state_solver.py
@@ -0,0 +1,152 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function, division
+
+# This program is part of 'MOOSE', the
+# Messaging Object Oriented Simulation Environment.
+#           Copyright (C) 2013 Upinder S. Bhalla. and NCBS
+# It is made available under the terms of the
+# GNU Lesser General Public License version 2.1
+# See the file COPYING.LIB for the full notice.
+# Monday 17 September 2018 01:49:30 PM IST
+# - Converted to a test script
+
+import math
+import numpy as np
+import moose
+print( "[INFO ] Using moose from %s" % moose.__file__ )
+
+def main():
+    compartment = makeModel()
+    ksolve = moose.Ksolve( '/model/compartment/ksolve' )
+    stoich = moose.Stoich( '/model/compartment/stoich' )
+    stoich.compartment = compartment
+    stoich.ksolve = ksolve
+    stoich.path = "/model/compartment/##"
+    state = moose.SteadyState( '/model/compartment/state' )
+
+    moose.reinit()
+    state.stoich = stoich
+    state.convergenceCriterion = 1e-6
+    moose.seed( 111 ) # Used when generating the samples in state space
+
+    b = moose.element( '/model/compartment/b' )
+    a = moose.element( '/model/compartment/a' )
+    c = moose.element( '/model/compartment/c' )
+    a.concInit = 0.1
+    deltaA = 0.002
+    num = 150
+    avec = []
+    bvec = []
+    moose.reinit()
+
+    # Now go up.
+    for i in range( 0, num ):
+        moose.start( 1.0 ) # Run the model for 1 seconds.
+        state.settle() # This function finds the steady states.
+        avec.append( a.conc )
+        bvec.append( b.conc )
+        a.concInit += deltaA
+
+    aa, bb = avec, bvec
+    got = np.mean(aa), np.std(aa)
+    expected = 0.24899, 0.08660
+    assert np.isclose(got, expected, atol = 1e-4).all(), "Got %s, expected %s" % (got, expected)
+    print( "[INFO ] Test 1 PASSED" )
+    
+
+    # Now go down.
+    avec = []
+    bvec = []
+    for i in range( 0, num ):
+        moose.start( 1.0 ) # Run the model for 1 seconds.
+        state.settle() # This function finds the steady states.
+        avec.append( a.conc )
+        bvec.append( b.conc )
+        a.concInit -= deltaA
+
+    aa, bb = avec,  bvec
+    got = np.mean(aa), np.std(aa)
+    expected = 0.251, 0.0866
+    assert np.isclose(got, expected, atol = 1e-4).all(), "Got %s, expected %s" % (got, expected)
+    print( "[INFO ] Test 2 PASSED" )
+
+    # Now aim for the middle. We do this by judiciously choosing a 
+    # start point that should be closer to the unstable fixed point.
+    avec = []
+    bvec = []
+    a.concInit = 0.28
+    b.conc = 0.15
+    for i in range( 0, 65 ):
+        moose.start( 1.0 ) # Run the model for 1 seconds.
+        state.settle() # This function finds the steady states.
+        avec.append( a.conc )
+        bvec.append( b.conc )
+        a.concInit -= deltaA
+
+    aa, bb = avec, bvec
+    got = np.mean(aa), np.std(aa)
+    expected = 0.216, 0.03752
+    assert np.isclose(got, expected, atol = 1e-4).all(), "Got %s, expected %s" % (got, expected)
+    print( "[INFO ] Test 3 PASSED" )
+    quit()
+
+def makeModel():
+    """ This function creates a bistable reaction system using explicit
+    MOOSE calls rather than load from a file.
+    The reaction is::
+
+        a ---b---> 2b    # b catalyzes a to form more of b.
+        2b ---c---> a    # c catalyzes b to form a.
+        a <======> 2b    # a interconverts to b.
+
+    """
+    # create container for model
+    model = moose.Neutral( 'model' )
+    compartment = moose.CubeMesh( '/model/compartment' )
+    compartment.volume = 1e-15
+    # the mesh is created automatically by the compartment
+    mesh = moose.element( '/model/compartment/mesh' ) 
+
+    # create molecules and reactions
+    a = moose.BufPool( '/model/compartment/a' )
+    b = moose.Pool( '/model/compartment/b' )
+    c = moose.Pool( '/model/compartment/c' )
+    enz1 = moose.Enz( '/model/compartment/b/enz1' )
+    enz2 = moose.Enz( '/model/compartment/c/enz2' )
+    cplx1 = moose.Pool( '/model/compartment/b/enz1/cplx' )
+    cplx2 = moose.Pool( '/model/compartment/c/enz2/cplx' )
+    reac = moose.Reac( '/model/compartment/reac' )
+
+    # connect them up for reactions
+    moose.connect( enz1, 'sub', a, 'reac' )
+    moose.connect( enz1, 'prd', b, 'reac' )
+    moose.connect( enz1, 'prd', b, 'reac' ) # Note 2 molecules of b.
+    moose.connect( enz1, 'enz', b, 'reac' )
+    moose.connect( enz1, 'cplx', cplx1, 'reac' )
+
+    moose.connect( enz2, 'sub', b, 'reac' )
+    moose.connect( enz2, 'sub', b, 'reac' ) # Note 2 molecules of b.
+    moose.connect( enz2, 'prd', a, 'reac' )
+    moose.connect( enz2, 'enz', c, 'reac' )
+    moose.connect( enz2, 'cplx', cplx2, 'reac' )
+
+    moose.connect( reac, 'sub', a, 'reac' )
+    moose.connect( reac, 'prd', b, 'reac' )
+    moose.connect( reac, 'prd', b, 'reac' ) # Note 2 order in b.
+
+    # Assign parameters
+    a.concInit = 1
+    b.concInit = 0
+    c.concInit = 0.01
+    enz1.kcat = 0.4
+    enz1.Km = 4
+    enz2.kcat = 0.6
+    enz2.Km = 0.01
+    reac.Kf = 0.001
+    reac.Kb = 0.01
+
+    return compartment
+
+# Run the 'main' if this script is executed standalone.
+if __name__ == '__main__':
+    main()
diff --git a/scripts/pull_subtree.sh b/scripts/pull_subtree.sh
index 71d70c29dacf513748eb3a5140bf4c7d433c3d9f..4836384b41cf659d537b3f2c357f6bd66419f410 100755
--- a/scripts/pull_subtree.sh
+++ b/scripts/pull_subtree.sh
@@ -24,11 +24,11 @@ set -x
 if [[ `pwd` == *"/moose" ]]; then
 
     git subtree pull --prefix moose-core \
-        https://github.com/BhallaLab/moose-core chhennapoda --squash
+        https://github.com/BhallaLab/moose-core master --squash
     git subtree pull --prefix moose-examples \
-        https://github.com/BhallaLab/moose-examples chhennapoda --squash 
+        https://github.com/BhallaLab/moose-examples master --squash 
     git subtree pull --prefix moose-gui \
-        https://github.com/BhallaLab/moose-gui chhennapoda --squash
+        https://github.com/BhallaLab/moose-gui master --squash
 
 else
     echo "Run this script from top-level git directory."