diff --git a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
index 2923f1ea00128a7d322aa9ed54eb6c9dd293b803..5c90e10f9dbe24c85d64bdb015ff634c26f6ee4b 100644
--- a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
+++ b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
@@ -92,8 +92,6 @@ FEMDeformableBodyModel::configure(const std::string& configFileName)
 bool
 FEMDeformableBodyModel::initialize()
 {
-    auto physicsMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(this->getModelGeometry());
-
     // prerequisite of for successfully initializing
     if (!m_forceModelGeometry || !m_forceModelConfiguration)
     {
@@ -101,16 +99,45 @@ FEMDeformableBodyModel::initialize()
         return false;
     }
 
+    auto physicsMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(this->getModelGeometry());
     m_vegaPhysicsMesh = physicsMesh->getAttachedVegaMesh();
 
-    this->initializeForceModel();
-    this->initializeMassMatrix();
-    this->initializeDampingMatrix();
-    this->initializeTangentStiffness();
+    if (!this->initializeForceModel())
+    {
+        return false;
+    }
+    ;
+    if (!this->initializeMassMatrix())
+    {
+        return false;
+    }
+    ;
+    if (!this->initializeDampingMatrix())
+    {
+        return false;
+    }
+    ;
+    if (!this->initializeTangentStiffness())
+    {
+        return false;
+    }
+    ;
     this->loadInitialStates();
-    this->loadBoundaryConditions();
-    this->initializeGravityForce();
-    this->initializeExplicitExternalForces();
+    if (!this->loadBoundaryConditions())
+    {
+        return false;
+    }
+    ;
+    if (!this->initializeGravityForce())
+    {
+        return false;
+    }
+    ;
+    if (!this->initializeExplicitExternalForces())
+    {
+        return false;
+    }
+    ;
 
     m_Feff.resize(m_numDOF);
     m_Finternal.resize(m_numDOF);
@@ -126,13 +153,18 @@ FEMDeformableBodyModel::initialize()
 void
 FEMDeformableBodyModel::loadInitialStates()
 {
+    if (m_numDOF == 0)
+    {
+        LOG(WARNING) << "FEMDeformableBodyModel::loadInitialStates() - Num. of degree of freedom is zero!";
+    }
+
     // For now the initial states are set to zero
     m_initialState = std::make_shared<kinematicState>(m_numDOF);
     m_previousState = std::make_shared<kinematicState>(m_numDOF);
     m_currentState = std::make_shared<kinematicState>(m_numDOF);
 }
 
-void
+bool
 FEMDeformableBodyModel::loadBoundaryConditions()
 {
     auto fileName = m_forceModelConfiguration->getStringOptionsMap().at("fixedDOFFilename");
@@ -140,7 +172,7 @@ FEMDeformableBodyModel::loadBoundaryConditions()
     if (fileName.empty())
     {
         LOG(WARNING) << "DeformableBodyModel::loadBoundaryConditions: The external boundary conditions file name is empty";
-        return;
+        return false;
     }
     else
     {
@@ -149,7 +181,6 @@ FEMDeformableBodyModel::loadBoundaryConditions()
         if (file.peek() == std::ifstream::traits_type::eof())
         {
             LOG(INFO) << "DeformableBodyModel::loadBoundaryConditions: The external boundary conditions file is empty";
-            return;
         }
 
         if (file.is_open())
@@ -167,6 +198,7 @@ FEMDeformableBodyModel::loadBoundaryConditions()
                 {
                     LOG(WARNING) << "FEMDeformableBodyModel::loadBoundaryConditions(): " <<
                         "The boundary condition node id provided is greater than number of nodes and hence excluded!!";
+                    return false;
                 }
             }
 
@@ -176,11 +208,13 @@ FEMDeformableBodyModel::loadBoundaryConditions()
         else
         {
             LOG(WARNING) << "DeformableBodyModel::loadBoundaryConditions: Could not open boundary conditions file!";
+            return false;
         }
     }
+    return true;
 }
 
-void
+bool
 FEMDeformableBodyModel::initializeForceModel()
 {
     const double g = m_forceModelConfiguration->getFloatsOptionsMap().at("gravity");
@@ -217,16 +251,19 @@ FEMDeformableBodyModel::initializeForceModel()
 
     default:
         LOG(WARNING) << "DeformableBodyModel::initializeForceModel: Unknown force model type";
+        return false;
     } //switch
+
+    return true;
 }
 
-void
+bool
 FEMDeformableBodyModel::initializeMassMatrix(const bool saveToDisk /*= false*/)
 {
     if (!m_forceModelGeometry)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeMassMatrix Force model geometry not set!";
-        return;
+        return false;
     }
 
     vega::SparseMatrix *vegaMatrix;
@@ -237,9 +274,11 @@ FEMDeformableBodyModel::initializeMassMatrix(const bool saveToDisk /*= false*/)
     this->initializeEigenMatrixFromVegaMatrix(*vegaMatrix, m_M);
 
     // TODO: Add option to save mass matrix to file
+
+    return true;
 }
 
-void
+bool
 FEMDeformableBodyModel::initializeDampingMatrix()
 {
     auto dampingLaplacianCoefficient = m_forceModelConfiguration->getFloatsOptionsMap().at("dampingLaplacianCoefficient");
@@ -250,13 +289,13 @@ FEMDeformableBodyModel::initializeDampingMatrix()
 
     if (!m_damped)
     {
-        return;
+        return false;
     }
 
-    if (dampingLaplacianCoefficient <= 0.0)
+    if (dampingLaplacianCoefficient < 0.0)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeDampingMatrix Damping coefficient is negative!";
-        return;
+        return false;
     }
 
     auto imstkVolMesh = std::static_pointer_cast<VolumetricMesh>(m_forceModelGeometry);
@@ -267,7 +306,7 @@ FEMDeformableBodyModel::initializeDampingMatrix()
     if (!meshGraph)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeDampingMatrix: Mesh graph not avaliable!";
-        return;
+        return false;
     }
 
     vega::SparseMatrix *matrix;
@@ -276,7 +315,7 @@ FEMDeformableBodyModel::initializeDampingMatrix()
     if (!matrix)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeDampingMatrix: Mesh Laplacian not avaliable!";
-        return;
+        return false;
     }
 
     matrix->ScalarMultiply(dampingLaplacianCoefficient);
@@ -284,15 +323,17 @@ FEMDeformableBodyModel::initializeDampingMatrix()
     m_vegaDampingMatrix.reset(matrix);
 
     this->initializeEigenMatrixFromVegaMatrix(*matrix, m_C);
+
+    return true;
 }
 
-void
+bool
 FEMDeformableBodyModel::initializeTangentStiffness()
 {
     if (!m_internalForceModel)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeTangentStiffness: Tangent stiffness cannot be initialized without force model";
-        return;
+        return false;
     }
 
     vega::SparseMatrix *matrix;
@@ -301,13 +342,13 @@ FEMDeformableBodyModel::initializeTangentStiffness()
     if (!matrix)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeTangentStiffness - Tangent stiffness matrix topology not avaliable!";
-        return;
+        return false;
     }
 
     if (!m_vegaMassMatrix)
     {
         LOG(WARNING) << "DeformableBodyModel::initializeTangentStiffness - Vega mass matrix doesn't exist!";
-        return;
+        return false;
     }
 
     matrix->BuildSubMatrixIndices(*m_vegaMassMatrix.get());
@@ -331,9 +372,11 @@ FEMDeformableBodyModel::initializeTangentStiffness()
     }
 
     m_internalForceModel->setTangentStiffness(m_vegaTangentStiffnessMatrix);
+
+    return true;
 }
 
-void
+bool
 FEMDeformableBodyModel::initializeGravityForce()
 {
     m_Fgravity.resize(m_numDOF);
@@ -341,6 +384,8 @@ FEMDeformableBodyModel::initializeGravityForce()
     const double gravity = m_forceModelConfiguration->getFloatsOptionsMap().at("gravity");
 
     m_vegaPhysicsMesh->computeGravity(m_Fgravity.data(), gravity);
+
+    return true;
 }
 
 void
@@ -450,11 +495,13 @@ FEMDeformableBodyModel::computeImplicitSystemLHS(const kinematicState& stateAtT,
     }
 }
 
-void
+bool
 FEMDeformableBodyModel::initializeExplicitExternalForces()
 {
     m_FexplicitExternal.resize(m_numDOF);
     m_FexplicitExternal.setZero();
+
+    return true;
 }
 
 void
diff --git a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
index 96d682189c6f1fd4ffd6853c855e832fc22316b5..b6ac2850ff78859c69095c066732724710fec8e3 100644
--- a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
+++ b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
@@ -110,37 +110,37 @@ public:
     ///
     /// \brief Load the boundary conditions from external file
     ///
-    void loadBoundaryConditions();
+    bool loadBoundaryConditions();
 
     ///
     /// \brief Initialize the force model
     ///
-    void initializeForceModel();
+    bool initializeForceModel();
 
     ///
     /// \brief Initialize the mass matrix from the mesh
     ///
-    void initializeMassMatrix(const bool saveToDisk = false);
+    bool initializeMassMatrix(const bool saveToDisk = false);
 
     ///
     /// \brief Initialize the damping (combines structural and viscous damping) matrix
     ///
-    void initializeDampingMatrix();
+    bool initializeDampingMatrix();
 
     ///
     /// \brief Initialize the tangent stiffness matrix
     ///
-    void initializeTangentStiffness();
+    bool initializeTangentStiffness();
 
     ///
     /// \brief Initialize the gravity force
     ///
-    void initializeGravityForce();
+    bool initializeGravityForce();
 
     ///
     /// \brief Initialize explicit external forces
     ///
-    void initializeExplicitExternalForces();
+    bool initializeExplicitExternalForces();
 
     ///
     /// \brief Initialize the Eigen matrix with data inside vega sparse matrix
diff --git a/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index 03c33944251207ed13b554327f8ddc2280a1c79b..258a189ebed14aab017bac164d6c58683eef665d 100644
--- a/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -147,22 +147,35 @@ PbdModel::initialize()
 
             if (strncmp("Corotation", &s[pos], len) == 0)
             {
-                LOG(INFO) << "Creating Corotation constraints";
-                this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::Corotation);
+                if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::Corotation))
+                {
+                    return false;
+                }
+                ;
             }
             else if (strncmp("NeoHookean", &s[pos], len) == 0)
             {
-                LOG(INFO) << "Creating Neohookean constraints";
-                this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::NeoHookean);
+                if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::NeoHookean))
+                {
+                    return false;
+                }
+                ;
             }
             else if (strncmp("Stvk", &s[pos], len) == 0)
             {
-                LOG(INFO) << "Creating StVenant-Kirchhoff constraints";
-                this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
+                if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK))
+                {
+                    return false;
+                }
+                ;
             }
             else
             { // default
-                this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
+                if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK))
+                {
+                    return false;
+                }
+                ;
             }
 
             float YoungModulus, PoissonRatio;
@@ -173,40 +186,56 @@ PbdModel::initialize()
         {
             float stiffness;
             sscanf(&s[len + 1], "%f", &stiffness);
-            LOG(INFO) << "Creating Volume constraints " << stiffness;
-            this->initializeVolumeConstraints(stiffness);
+            if (!this->initializeVolumeConstraints(stiffness))
+            {
+                return false;
+            }
+            ;
         }
         else if (strncmp("Distance", &s[0], len) == 0)
         {
             float stiffness;
             sscanf(&s[len + 1], "%f", &stiffness);
-            LOG(INFO) << "Creating Distance constraints " << stiffness;
-            this->initializeDistanceConstraints(stiffness);
+            if (!this->initializeDistanceConstraints(stiffness))
+            {
+                return false;
+            }
+            ;
         }
         else if (strncmp("Area", &s[0], len) == 0)
         {
             float stiffness;
             sscanf(&s[len + 1], "%f", &stiffness);
-            LOG(INFO) << "Creating Area constraints " << stiffness;
-            this->initializeAreaConstraints(stiffness);
+            if (!this->initializeAreaConstraints(stiffness))
+            {
+                return false;
+            }
+            ;
         }
         else if (strncmp("Dihedral", &s[0], len) == 0)
         {
             float stiffness;
             sscanf(&s[len + 1], "%f", &stiffness);
-            LOG(INFO) << "Creating Dihedral constraints " << stiffness;
-            this->initializeDihedralConstraints(stiffness);
+            if (!this->initializeDihedralConstraints(stiffness))
+            {
+                return false;
+            }
+            ;
         }
         else if (strncmp("ConstantDensity", &s[0], len) == 0)
         {
             float stiffness;
             sscanf(&s[len + 1], "%f", &stiffness);
-            LOG(INFO) << "Creating Constant Density constraints ";
-            this->initializeConstantDensityConstraint(stiffness);
+            if (!this->initializeConstantDensityConstraint(stiffness))
+            {
+                return false;
+            }
+            ;
         }
         else
         {
-            exit(0);
+            LOG(WARNING) << "PbdModel::initialize() - Type of PBD constraints not identified!";
+            return false;
         }
     }
     return true;
diff --git a/Base/Scene/imstkScene.cpp b/Base/Scene/imstkScene.cpp
index 7585394fe8626007ac35c09107078f296459cd70..2a3377ef79474123de880682ed31c12f849ab62d 100644
--- a/Base/Scene/imstkScene.cpp
+++ b/Base/Scene/imstkScene.cpp
@@ -27,6 +27,22 @@
 
 namespace imstk
 {
+
+bool
+Scene::initialize()
+{
+    for (auto const& it : m_sceneObjectsMap)
+    {
+        auto sceneObject = it.second;
+        if (!sceneObject->initialize())
+        {
+            LOG(WARNING) << "Error initializing scene object: " << sceneObject->getName();
+            return false;
+        }
+    }
+    return true;
+}
+
 bool
 Scene::isObjectRegistered(std::string sceneObjectName) const
 {
diff --git a/Base/Scene/imstkScene.h b/Base/Scene/imstkScene.h
index 426d8e9783a28a7ef72f1b3428ee06af7ce60149..4443c8bdc49a8b4a3a7e450b9a0a35ffa45341b8 100644
--- a/Base/Scene/imstkScene.h
+++ b/Base/Scene/imstkScene.h
@@ -57,6 +57,11 @@ public:
     ///
     ~Scene() = default;
 
+    ///
+    /// \brief Initialize the scene
+    ///
+    bool initialize();
+
     ///
     /// \brief Returns true if the object with a given name is registered, else false
     ///
diff --git a/Base/SceneElements/Objects/imstkCollidingObject.h b/Base/SceneElements/Objects/imstkCollidingObject.h
index 5a1fed3aafab53769eaf12349bff762edbbae256..1a4c343c45f5bf4ce649eed6e236bbf573ac4322 100644
--- a/Base/SceneElements/Objects/imstkCollidingObject.h
+++ b/Base/SceneElements/Objects/imstkCollidingObject.h
@@ -78,6 +78,21 @@ public:
     ///
     virtual void updateGeometries() override;
 
+    ///
+    /// \brief Initialize the scene object
+    ///
+    virtual bool initialize()
+    {
+        if (SceneObject::initialize())
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
 protected:
 
     std::shared_ptr<Geometry> m_collidingGeometry;       ///> Geometry for collisions
diff --git a/Base/SceneElements/Objects/imstkDeformableObject.cpp b/Base/SceneElements/Objects/imstkDeformableObject.cpp
index 0e8ca53ead70d55d34e1261edd3ba8574c7c05cb..0b2ad8125b80e9d9a40178a21d34a6a031fea1ac 100644
--- a/Base/SceneElements/Objects/imstkDeformableObject.cpp
+++ b/Base/SceneElements/Objects/imstkDeformableObject.cpp
@@ -26,11 +26,9 @@ namespace imstk
 Vectord&
 DeformableObject::getContactForce()
 {
-    //m_defModel = std::dynamic_pointer_cast<imstk::DeformableBodyModel>(m_dynamicalModel);
-
     if (!m_defModel)
     {
-        LOG(WARNING) << "Dynamics pointer cast failure in DeformableObject::getContactForce()";
+        LOG(WARNING) << "deformation model pointer not valid DeformableObject::getContactForce()";
     }
 
     return m_defModel->getContactForce();
@@ -40,12 +38,15 @@ bool
 DeformableObject::initialize()
 {
     m_defModel = std::dynamic_pointer_cast<FEMDeformableBodyModel>(m_dynamicalModel);
-    if (m_defModel == nullptr)
+    if (m_defModel)
+    {
+        return DynamicObject::initialize();
+    }
+    else
     {
-        LOG(WARNING) << "Dynamic model set is not of expected type (DeformableBodyModel)!";
+        LOG(WARNING) << "Dynamics pointer cast failure in DeformableObject::initialize()";
         return false;
     }
-    return true;
 }
 
 const Vectord&
diff --git a/Base/SceneElements/Objects/imstkDeformableObject.h b/Base/SceneElements/Objects/imstkDeformableObject.h
index 24d3468043c46e39c17ce5dd976a1701ca28fb1f..4797896c6533b9f4ae45d4225d1b8efce7703b13 100644
--- a/Base/SceneElements/Objects/imstkDeformableObject.h
+++ b/Base/SceneElements/Objects/imstkDeformableObject.h
@@ -56,9 +56,9 @@ public:
     ~DeformableObject() = default;
 
     ///
-    /// \brief
+    /// \brief Initialize the deformable object
     ///
-    bool initialize();
+    bool initialize() override;
 
     ///
     /// \brief Initialize the kinematic state of the body
diff --git a/Base/SceneElements/Objects/imstkDynamicObject.h b/Base/SceneElements/Objects/imstkDynamicObject.h
index 6f6aeb561ea28bb0315e0912ae20d4d636514440..e6010bde8dae62c5797beb5466d1f604d41e3863 100644
--- a/Base/SceneElements/Objects/imstkDynamicObject.h
+++ b/Base/SceneElements/Objects/imstkDynamicObject.h
@@ -114,6 +114,21 @@ public:
         }
     }
 
+    ///
+    /// \brief Initialize the scene object
+    ///
+    virtual bool initialize() override
+    {
+        if (CollidingObject::initialize())
+        {
+            return m_dynamicalModel->initialize();
+        }
+        else
+        {
+            return false;
+        }
+    }
+
 protected:
 
     ///
@@ -122,12 +137,12 @@ protected:
     DynamicObject(std::string name) : CollidingObject(name){}
 
     std::shared_ptr<DynamicalModel<StateType>> m_dynamicalModel;        ///> Dynamical model
-    std::shared_ptr<Geometry> m_physicsGeometry;                ///> Geometry used for Physics
+    std::shared_ptr<Geometry> m_physicsGeometry;                        ///> Geometry used for Physics
 
     //Maps
-    std::shared_ptr<GeometryMap> m_physicsToCollidingGeomMap;   ///> Maps from Physics to collision geometry
-    std::shared_ptr<GeometryMap> m_physicsToVisualGeomMap;      ///> Maps from Physics to visual geometry
-    bool m_updateVisualFromPhysicsGeometry = true; ///> Defines if visual is updated from colliding mapping or physics mapping
+    std::shared_ptr<GeometryMap> m_physicsToCollidingGeomMap;           ///> Maps from Physics to collision geometry
+    std::shared_ptr<GeometryMap> m_physicsToVisualGeomMap;              ///> Maps from Physics to visual geometry
+    bool m_updateVisualFromPhysicsGeometry = true;                      ///> Defines if visual is updated from colliding mapping or physics mapping
 };
 } // imstk
 
diff --git a/Base/SceneElements/Objects/imstkPbdObject.cpp b/Base/SceneElements/Objects/imstkPbdObject.cpp
index fbbdc70d0273ca533b71c69ca9b492711d0b6420..61537799d503cae7a14e12f517073a48aa2f0b33 100644
--- a/Base/SceneElements/Objects/imstkPbdObject.cpp
+++ b/Base/SceneElements/Objects/imstkPbdObject.cpp
@@ -27,155 +27,18 @@
 namespace imstk
 {
 bool
-PbdObject::initialize(int nCons, ...)
+PbdObject::initialize()
 {
-    std::shared_ptr<PointSet> mesh = std::dynamic_pointer_cast<PointSet>(m_physicsGeometry);
-    if (mesh == nullptr)
-    {
-        LOG(WARNING) << "Physics geometry is not a mesh!";
-        return false;
-    }
-
     m_pbdModel = std::dynamic_pointer_cast<PbdModel>(m_dynamicalModel);
-    if (m_pbdModel == nullptr)
-    {
-        LOG(WARNING) << "Dynamic model set is not of expected type (PbdModel)!";
-        return false;
-    }
-
-    m_pbdModel->setModelGeometry(mesh);
-    m_pbdModel->initialize();
-
-    va_list args;
-    va_start(args, nCons);
-    for(int i = 0; i < nCons; ++i)
-    {
-        char* s = va_arg(args,char*);
-        int len = 0;
-        while (s[len] != ' ' && s[len] != '\0')
-        {
-            ++len;
-        }
-
-        if (strncmp("FEM",&s[0],len) == 0)
-        {
-            int pos = len+1;
-            len = 0;
-            while (s[pos+len] != ' ' && s[pos+len] != '\0')
-            {
-                ++len;
-            }
-
-            if (strncmp("Corotation",&s[pos],len) == 0)
-            {
-                LOG(INFO) << "Creating Corotation constraints";
-                m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::Corotation);
-            }
-            else if (strncmp("NeoHookean",&s[pos],len) == 0)
-            {
-                LOG(INFO) << "Creating Neohookean constraints";
-                m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::NeoHookean);
-            }
-            else if (strncmp("Stvk",&s[pos],len) == 0)
-            {
-                LOG(INFO) << "Creating StVenant-Kirchhoff constraints";
-                m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
-            }
-            else
-            { // default
-                m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
-            }
-
-            float YoungModulus, PoissonRatio;
-            sscanf(&s[pos+len+1], "%f %f", &YoungModulus, &PoissonRatio);
-            m_pbdModel->computeLameConstants(YoungModulus, PoissonRatio);
-        }
-        else if (strncmp("Volume",&s[0],len) == 0)
-        {
-            float stiffness;
-            sscanf(&s[len+1], "%f", &stiffness);
-            LOG(INFO) << "Creating Volume constraints " << stiffness;
-            m_pbdModel->initializeVolumeConstraints(stiffness);
-        }
-        else if (strncmp("Distance",&s[0],len) == 0)
-        {
-            float stiffness;
-            sscanf(&s[len+1], "%f", &stiffness);
-            LOG(INFO) << "Creating Distance constraints " << stiffness;
-            m_pbdModel->initializeDistanceConstraints(stiffness);
-        }
-        else if (strncmp("Area",&s[0],len) == 0)
-        {
-            float stiffness;
-            sscanf(&s[len+1], "%f", &stiffness);
-            LOG(INFO) << "Creating Area constraints " << stiffness;
-            m_pbdModel->initializeAreaConstraints(stiffness);
-        }
-        else if (strncmp("Dihedral",&s[0],len) == 0)
-        {
-            float stiffness;
-            sscanf(&s[len+1], "%f", &stiffness);
-            LOG(INFO) << "Creating Dihedral constraints " << stiffness;
-            m_pbdModel->initializeDihedralConstraints(stiffness);
-        }
-        else if (strncmp("ConstantDensity", &s[0], len) == 0)
-        {
-            float stiffness;
-            sscanf(&s[len + 1], "%f", &stiffness);
-            LOG(INFO) << "Creating Constant Density constraints ";
-            m_pbdModel->initializeConstantDensityConstraint(stiffness);
-        }
-        else
-        {
-            exit(0);
-        }
-    }
-
-    m_pbdModel->setUniformMass(va_arg(args, double));
-
-    if (nCons > 0)
+    if (m_pbdModel)
     {
-        char* gstring = va_arg(args,char*);
-        float x,y,z;
-        sscanf(gstring,"%f %f %f", &x, &y, &z);
-        Vec3d g(x,y,z);
-        m_pbdModel->setGravity(g);
-        m_pbdModel->setTimeStep(va_arg(args, double));
-        char *s = va_arg(args,char*);
-
-        if (strlen(s) > 0)
-        {
-            while (1)
-            {
-                int idx = atoi(s);
-                m_pbdModel->setFixedPoint(idx - 1);
-                while (*s != ' ' && *s != '\0')
-                {
-                    ++s;
-                }
-
-                if (*s == '\0')
-                {
-                    break;
-                }
-                else
-                {
-                    ++s;
-                }
-            }
-        }
-        m_pbdModel->setMaxNumIterations(va_arg(args,int));
+        return DynamicObject::initialize();
     }
-
-    if (m_collidingGeometry)
+    else
     {
-        m_pbdModel->setProximity(va_arg(args,double));
-        m_pbdModel->setContactStiffness(va_arg(args,double));
+        LOG(WARNING) << "Dynamics pointer cast failure in PbdObject::initialize()";
+        return false;
     }
-
-    m_pbdModel->setNumDegreeOfFreedom(m_pbdModel->getModelGeometry()->getNumVertices() * 3);
-
-    return true;
 }
 
 void
diff --git a/Base/SceneElements/Objects/imstkPbdObject.h b/Base/SceneElements/Objects/imstkPbdObject.h
index b8f06e0d2f76a3771e58e8b98ba5434e2ad6ed1c..554f3fd6bba673ad70fe19c44a68ef0723182be3 100644
--- a/Base/SceneElements/Objects/imstkPbdObject.h
+++ b/Base/SceneElements/Objects/imstkPbdObject.h
@@ -56,10 +56,9 @@ public:
     virtual ~PbdObject() = default;
 
     ///
-    /// \brief Initialize the pbd configuration
-    /// TODO: Parse from config file
+    /// \brief Initialize the pbd scene object
     ///
-    bool initialize(int nCons, ...);
+    bool initialize() override;
 
     ///
     /// \brief Update the position based on Verlet time stepping rule
diff --git a/Base/SceneElements/Objects/imstkSceneObject.h b/Base/SceneElements/Objects/imstkSceneObject.h
index 6b9655d007bd4fb445cacf584d44f4f2e17086ce..4e9b24a1967e390fea28ce79c295770c1ef5ef66 100644
--- a/Base/SceneElements/Objects/imstkSceneObject.h
+++ b/Base/SceneElements/Objects/imstkSceneObject.h
@@ -88,6 +88,11 @@ public:
     ///
     virtual void updateGeometries() {};
 
+    ///
+    /// \brief Initialize the scene object
+    ///
+    virtual bool initialize(){ return true; };
+
 protected:
     ///
     /// \brief Assigns the type of the object