diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp
index c897a9e6cf..78980a9a71 100644
--- a/nestkernel/nest_names.cpp
+++ b/nestkernel/nest_names.cpp
@@ -488,6 +488,8 @@ const Name stimulus_source( "stimulus_source" );
 const Name stop( "stop" );
 const Name structural_plasticity_synapses( "structural_plasticity_synapses" );
 const Name structural_plasticity_update_interval( "structural_plasticity_update_interval" );
+const Name structural_plasticity_gaussian_kernel_sigma( "structural_plasticity_gaussian_kernel_sigma" );
+const Name structural_plasticity_cache_probabilities( "structural_plasticity_cache_probabilities" );
 const Name surrogate_gradient( "surrogate_gradient" );
 const Name surrogate_gradient_function( "surrogate_gradient_function" );
 const Name synapse_id( "synapse_id" );
diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h
index b6d964e04c..fa991f7b23 100644
--- a/nestkernel/nest_names.h
+++ b/nestkernel/nest_names.h
@@ -516,6 +516,8 @@ extern const Name stimulus_source;
 extern const Name stop;
 extern const Name structural_plasticity_synapses;
 extern const Name structural_plasticity_update_interval;
+extern const Name structural_plasticity_gaussian_kernel_sigma;
+extern const Name structural_plasticity_cache_probabilities;
 extern const Name surrogate_gradient;
 extern const Name surrogate_gradient_function;
 extern const Name synapse_id;
diff --git a/nestkernel/sp_manager.cpp b/nestkernel/sp_manager.cpp
index 4705b7de64..83785e8718 100644
--- a/nestkernel/sp_manager.cpp
+++ b/nestkernel/sp_manager.cpp
@@ -24,6 +24,7 @@
 
 // C++ includes:
 #include <algorithm>
+#include <random>
 
 // Includes from nestkernel:
 #include "conn_builder.h"
@@ -33,6 +34,7 @@
 #include "kernel_manager.h"
 #include "nest_names.h"
 #include "sp_manager_impl.h"
+#include "spatial.h"
 
 namespace nest
 {
@@ -41,6 +43,8 @@ SPManager::SPManager()
   : ManagerInterface()
   , structural_plasticity_update_interval_( 10000. )
   , structural_plasticity_enabled_( false )
+  , structural_plasticity_gaussian_kernel_sigma_( -1 )
+  , structural_plasticity_cache_probabilities_( false )
   , sp_conn_builders_()
   , growthcurve_factories_()
   , growthcurvedict_( new Dictionary() )
@@ -64,6 +68,8 @@ SPManager::initialize( const bool adjust_number_of_threads_or_rng_only )
 
   structural_plasticity_update_interval_ = 10000.;
   structural_plasticity_enabled_ = false;
+  structural_plasticity_gaussian_kernel_sigma_ = -1;
+  structural_plasticity_cache_probabilities_ = false;
 }
 
 void
@@ -106,6 +112,9 @@ SPManager::get_status( DictionaryDatum& d )
   }
 
   def< double >( d, names::structural_plasticity_update_interval, structural_plasticity_update_interval_ );
+  def< double >( d, names::structural_plasticity_gaussian_kernel_sigma, structural_plasticity_gaussian_kernel_sigma_ );
+  def< bool >( d, names::structural_plasticity_cache_probabilities, structural_plasticity_cache_probabilities_ );
+
 
   ArrayDatum growth_curves;
   for ( auto const& element : *growthcurvedict_ )
@@ -119,6 +128,10 @@ void
 SPManager::set_status( const DictionaryDatum& d )
 {
   updateValue< double >( d, names::structural_plasticity_update_interval, structural_plasticity_update_interval_ );
+  updateValue< double >(
+    d, names::structural_plasticity_gaussian_kernel_sigma, structural_plasticity_gaussian_kernel_sigma_ );
+  updateValue< bool >(
+    d, names::structural_plasticity_cache_probabilities, structural_plasticity_cache_probabilities_ );
 
   if ( not d->known( names::structural_plasticity_synapses ) )
   {
@@ -168,6 +181,210 @@ SPManager::set_status( const DictionaryDatum& d )
   }
 }
 
+void
+SPManager::gather_global_positions_and_ids()
+{
+  std::vector< double > local_positions;
+  std::vector< int > local_ids;
+  std::vector< int > displacements;
+
+  // Collect local positions and IDs
+  for ( size_t tid = 0; tid < kernel().vp_manager.get_num_threads(); ++tid )
+  {
+    const SparseNodeArray& local_nodes = kernel().node_manager.get_local_nodes( tid );
+
+    for ( auto node_it = local_nodes.begin(); node_it < local_nodes.end(); ++node_it )
+    {
+      int node_id = node_it->get_node_id();
+      if ( node_id < 1 )
+      {
+        throw std::runtime_error( "Invalid neuron ID (must be >= 1)." );
+      }
+
+      std::vector< double > pos = get_position( node_id );
+
+      if ( std::none_of( pos.begin(), pos.end(), []( double v ) { return std::isnan( v ); } ) )
+      {
+        local_ids.push_back( node_id );
+        local_positions.insert( local_positions.end(), pos.begin(), pos.end() );
+      }
+    }
+  }
+
+  // Communicate positions and IDs
+  kernel().mpi_manager.communicate( local_positions, global_positions, displacements );
+  kernel().mpi_manager.communicate( local_ids, global_ids, displacements );
+
+  // Validate global_positions size consistency with global_ids
+  size_t num_neurons = global_ids.size();
+  size_t total_positions = global_positions.size();
+
+  if ( num_neurons == 0 )
+  {
+    throw std::runtime_error( "No neurons found. Please provide positions, or disable distance dependency" );
+  }
+  if ( total_positions == 0 )
+  {
+    throw std::runtime_error( "No positions found. Please provide positions, or disable distance dependency." );
+  }
+
+  if ( total_positions % num_neurons != 0 )
+  {
+    throw std::runtime_error( "Mismatch in global positions dimensionality." );
+  }
+
+  pos_dim = total_positions / num_neurons;
+
+  // Pair global_ids with their positions
+  std::vector< std::pair< int, std::vector< double > > > id_pos_pairs;
+  id_pos_pairs.reserve( num_neurons );
+  for ( size_t i = 0; i < num_neurons; ++i )
+  {
+    int node_id = global_ids[ i ];
+    std::vector< double > pos( global_positions.begin() + i * pos_dim, global_positions.begin() + ( i + 1 ) * pos_dim );
+    id_pos_pairs.emplace_back( node_id, pos );
+  }
+
+  // Sort id_pos_pairs based on node_id to ensure ordering from 1 to num_neurons
+  std::sort( id_pos_pairs.begin(),
+    id_pos_pairs.end(),
+    []( const std::pair< int, std::vector< double > >& a, const std::pair< int, std::vector< double > >& b ) -> bool
+    { return a.first < b.first; } );
+
+  // Verify that IDs are sequential
+  for ( size_t i = 0; i < num_neurons; ++i )
+  {
+    if ( id_pos_pairs[ i ].first != static_cast< int >( i + 1 ) )
+    {
+      throw std::runtime_error( "Neuron IDs are not sequential after sorting." );
+    }
+  }
+
+  // Assign sorted positions to temp_positions
+  std::vector< double > temp_positions( num_neurons * pos_dim, 0.0 );
+  for ( size_t i = 0; i < num_neurons; ++i )
+  {
+    std::copy( id_pos_pairs[ i ].second.begin(), id_pos_pairs[ i ].second.end(), temp_positions.begin() + i * pos_dim );
+  }
+
+  // Update global_positions with sorted positions
+  global_positions = std::move( temp_positions );
+}
+
+// This method uses a formula based on triangular numbers
+// to map two ids to one index indepndent of theirs order
+int
+SPManager::get_neuron_pair_index( int id1, int id2 )
+{
+  int max_id = std::max( id1, id2 );
+  int min_id = std::min( id1, id2 );
+  int index = ( ( max_id ) * ( max_id - 1 ) ) / 2 + ( min_id - 1 );
+  return index;
+}
+
+
+// Method to perform roulette wheel selection
+int
+SPManager::roulette_wheel_selection( const std::vector< double >& probabilities, double rnd )
+{
+  if ( probabilities.empty() )
+  {
+    throw std::runtime_error( "Probabilities vector is empty." );
+  }
+
+  std::vector< double > cumulative( probabilities.size() );
+  std::partial_sum( probabilities.begin(), probabilities.end(), cumulative.begin() );
+
+  // Ensure the sum of probabilities is greater than zero
+  double sum = cumulative.back();
+  if ( sum < 0.0 )
+  {
+    throw std::runtime_error( "Sum of probabilities must be greater than zero." );
+  }
+
+
+  // Generate a random number in the range [0, sum)
+  double randomValue = rnd * sum;
+
+  // Perform binary search to find the selected index
+  auto it = std::lower_bound( cumulative.begin(), cumulative.end(), randomValue );
+  return static_cast< int >( std::distance( cumulative.begin(), it ) );
+}
+
+
+double
+SPManager::gaussian_kernel( const std::vector< double >& pos1, const std::vector< double >& pos2, const double sigma )
+{
+  double distanceSquared = 0.0;
+  for ( size_t i = 0; i < pos1.size(); ++i )
+  {
+    double diff = pos2[ i ] - pos1[ i ];
+    distanceSquared += diff * diff;
+  }
+  return std::exp( -distanceSquared / ( sigma * sigma ) );
+}
+
+void
+SPManager::build_probability_list()
+{
+  size_t num_neurons = global_ids.size();
+
+  if ( global_positions.size() % num_neurons != 0 )
+  {
+    throw std::runtime_error( "Mismatch in global positions dimensionality." );
+  }
+
+  // Resize the probability list to accommodate all neuron pairs.
+  size_t total_pairs = ( num_neurons * ( num_neurons + 1 ) ) / 2;
+  probability_list.resize( total_pairs, -1.0 );
+
+  // Calculate probabilities for connections between all pairs of neurons.
+  for ( size_t i = 0; i < num_neurons; ++i )
+  {
+    size_t id_i = i + 1;
+    if ( id_i < 1 || id_i > num_neurons )
+    {
+      std::cerr << "Error: Neuron ID " << id_i << " out of valid range." << std::endl;
+      continue;
+    }
+
+    std::vector< double > pos_i(
+      global_positions.begin() + pos_dim * ( id_i - 1 ), global_positions.begin() + pos_dim * id_i );
+
+    for ( size_t j = i; j < num_neurons; ++j )
+    {
+      size_t id_j = j + 1;
+      if ( id_j < 1 || id_j > num_neurons )
+      {
+        std::cerr << "Error: Neuron ID " << id_j << " out of valid range." << std::endl;
+        continue;
+      }
+
+      size_t index = get_neuron_pair_index( id_i, id_j );
+
+      if ( index >= probability_list.size() )
+      {
+        std::cerr << "Error: Index out of bounds: " << index << " for ids " << id_i << " and " << id_j << std::endl;
+        continue;
+      }
+
+      if ( id_i == id_j )
+      {
+        probability_list[ index ] = 0.0; // Assign zero probability for self-connections
+      }
+      else
+      {
+        std::vector< double > pos_j(
+          global_positions.begin() + pos_dim * ( id_j - 1 ), global_positions.begin() + pos_dim * id_j );
+
+        double prob = gaussian_kernel( pos_i, pos_j, structural_plasticity_gaussian_kernel_sigma_ );
+        probability_list[ index ] = prob;
+      }
+    }
+  }
+}
+
+
 long
 SPManager::builder_min_delay() const
 {
@@ -409,26 +626,39 @@ SPManager::create_synapses( std::vector< size_t >& pre_id,
   serialize_id( pre_id, pre_n, pre_id_rnd );
   serialize_id( post_id, post_n, post_id_rnd );
 
-  // Shuffle only the largest vector
-  if ( pre_id_rnd.size() > post_id_rnd.size() )
+  std::vector< size_t > pre_ids_results;
+  std::vector< size_t > post_ids_results;
+
+  if ( structural_plasticity_gaussian_kernel_sigma_ <= 0 )
   {
-    // we only shuffle the n first items,
-    // where n is the number of postsynaptic elements
-    global_shuffle( pre_id_rnd, post_id_rnd.size() );
-    pre_id_rnd.resize( post_id_rnd.size() );
+    // Shuffle only the largest vector
+    if ( pre_id_rnd.size() > post_id_rnd.size() )
+    {
+      // we only shuffle the n first items,
+      // where n is the number of postsynaptic elements
+      global_shuffle( pre_id_rnd, post_id_rnd.size() );
+      pre_id_rnd.resize( post_id_rnd.size() );
+    }
+    else
+    {
+      // we only shuffle the n first items,
+      // where n is the number of pre synaptic elements
+      global_shuffle( post_id_rnd, pre_id_rnd.size() );
+      post_id_rnd.resize( pre_id_rnd.size() );
+    }
+
+    pre_ids_results = pre_id_rnd;
+    post_ids_results = post_id_rnd;
   }
   else
   {
-    // we only shuffle the n first items,
-    // where n is the number of pre synaptic elements
-    global_shuffle( post_id_rnd, pre_id_rnd.size() );
-    post_id_rnd.resize( pre_id_rnd.size() );
+    global_shuffle_spatial( pre_id_rnd, post_id_rnd, pre_ids_results, post_ids_results );
   }
 
   // create synapse
-  sp_conn_builder->sp_connect( pre_id_rnd, post_id_rnd );
+  sp_conn_builder->sp_connect( pre_ids_results, post_ids_results );
 
-  return not pre_id_rnd.empty();
+  return not pre_ids_results.empty();
 }
 
 void
@@ -663,6 +893,87 @@ nest::SPManager::global_shuffle( std::vector< size_t >& v, size_t n )
   }
   v = v2;
 }
+void
+SPManager::global_shuffle_spatial( std::vector< size_t >& pre_ids,
+  std::vector< size_t >& post_ids,
+  std::vector< size_t >& pre_ids_results,
+  std::vector< size_t >& post_ids_results )
+{
+  size_t maxIterations = std::min( pre_ids.size(), post_ids.size() );
+
+  for ( size_t iteration = 0; iteration < maxIterations; ++iteration )
+  {
+    if ( pre_ids.empty() || post_ids.empty() )
+    {
+      break; // Stop if either vector is empty
+    }
+
+    size_t pre_id = pre_ids.back();
+    pre_ids.pop_back();
+
+    std::vector< double > probabilities;
+    std::vector< size_t > valid_post_ids;
+    double rnd;
+    for ( size_t post_id : post_ids )
+    {
+      if ( post_id == pre_id )
+      {
+        continue; // Skip self-connections
+      }
+
+      double prob;
+      if ( structural_plasticity_cache_probabilities_ )
+      {
+        // Retrieve cached probability for the neuron pair
+        int pair_index = get_neuron_pair_index( pre_id, post_id );
+        if ( pair_index < 0 || pair_index >= static_cast< int >( probability_list.size() ) )
+        {
+          std::cerr << "Error: index out of bounds for pair (" << pre_id << ", " << post_id << ")" << std::endl;
+          continue;
+        }
+        prob = probability_list[ pair_index ];
+      }
+      else
+      {
+        size_t pre_index = pre_id - 1;
+        std::vector< double > pre_pos(
+          global_positions.begin() + pre_index * pos_dim, global_positions.begin() + ( pre_index + 1 ) * pos_dim );
+
+        size_t post_index = post_id - 1;
+        std::vector< double > post_pos(
+          global_positions.begin() + post_index * pos_dim, global_positions.begin() + ( post_index + 1 ) * pos_dim );
+
+        prob = gaussian_kernel( pre_pos, post_pos, structural_plasticity_gaussian_kernel_sigma_ );
+      }
+      if ( prob > 0 )
+      {
+        probabilities.push_back( prob );
+        valid_post_ids.push_back( post_id );
+      }
+    }
+
+    if ( probabilities.empty() )
+    {
+      continue; // Skip if no valid connections are found
+    }
+
+    rnd = get_rank_synced_rng()->drand();
+
+    // Select a post-synaptic neuron using roulette wheel selection
+    int selected_post_idx = roulette_wheel_selection( probabilities, rnd );
+    size_t selected_post_id = valid_post_ids[ selected_post_idx ];
+
+    // Remove the selected post-synaptic neuron from the list
+    auto post_it = std::find( post_ids.begin(), post_ids.end(), selected_post_id );
+    if ( post_it != post_ids.end() )
+    {
+      post_ids.erase( post_it );
+    }
+
+    pre_ids_results.push_back( pre_id );
+    post_ids_results.push_back( selected_post_id );
+  }
+}
 
 
 void
@@ -685,6 +996,14 @@ nest::SPManager::enable_structural_plasticity()
       "has been set to false." );
   }
   structural_plasticity_enabled_ = true;
+  if ( structural_plasticity_gaussian_kernel_sigma_ > 0 )
+  {
+    gather_global_positions_and_ids();
+    if ( structural_plasticity_cache_probabilities_ )
+    {
+      build_probability_list();
+    }
+  }
 }
 
 void
@@ -693,4 +1012,4 @@ nest::SPManager::disable_structural_plasticity()
   structural_plasticity_enabled_ = false;
 }
 
-} // namespace nest
+} // namespace nest
\ No newline at end of file
diff --git a/nestkernel/sp_manager.h b/nestkernel/sp_manager.h
index 79cc92a72a..d2d8e66f1a 100644
--- a/nestkernel/sp_manager.h
+++ b/nestkernel/sp_manager.h
@@ -24,6 +24,7 @@
 #define SP_MANAGER_H
 
 // C++ includes:
+#include <random>
 #include <vector>
 
 // Includes from libnestutil:
@@ -137,6 +138,8 @@ class SPManager : public ManagerInterface
 
   double get_structural_plasticity_update_interval() const;
 
+  double get_structural_plasticity_gaussian_kernel_sigma() const;
+
   /**
    * Returns the minimum delay of all SP builders.
    *
@@ -187,6 +190,80 @@ class SPManager : public ManagerInterface
   void serialize_id( std::vector< size_t >& id, std::vector< int >& n, std::vector< size_t >& res );
   void global_shuffle( std::vector< size_t >& v );
   void global_shuffle( std::vector< size_t >& v, size_t n );
+  /**
+   * Calculate a unique index for a pair of neuron IDs for efficient lookup.
+   * The index is calculated independent of the order
+   *
+   * @param id1 First neuron ID.
+   * @param id2 Second neuron ID.
+   * @return Unique index corresponding to the neuron pair.
+   */
+  int get_neuron_pair_index( int id1, int id2 );
+
+  /**
+   * Compute the Gaussian kernel value between two positions.
+   *
+   * @param pos1 Position of the first neuron.
+   * @param pos2 Position of the second neuron.
+   * @param sigma Standard deviation for the Gaussian kernel.
+   * @return Gaussian kernel value.
+   */
+  double gaussian_kernel( const std::vector< double >& pos1, const std::vector< double >& pos2, const double sigma );
+
+  /**
+   * Perform global shuffling of pre- and post-synaptic neurons based on spatial probabilities.
+   *
+   * @param pre_ids Vector of pre-synaptic neuron IDs.
+   * @param post_ids Vector of post-synaptic neuron IDs.
+   * @param pre_ids_results Vector to store shuffled pre-synaptic IDs.
+   * @param post_ids_results Vector to store shuffled post-synaptic IDs.
+   */
+  void global_shuffle_spatial( std::vector< size_t >& pre_ids,
+    std::vector< size_t >& post_ids,
+    std::vector< size_t >& pre_ids_results,
+    std::vector< size_t >& post_ids_results );
+
+  /**
+   * Build a probability list for neuron connections based on spatial properties.
+   */
+  void build_probability_list();
+
+  /**
+   * Gather global neuron positions and IDs from all nodes.
+   */
+  void gather_global_positions_and_ids();
+
+  /**
+   * Perform roulette wheel selection to randomly select an index based on probabilities.
+   *
+   * @param probabilities Vector of probabilities for selection.
+   * @param rnd Random number.
+   * @return Selected index.
+   */
+  int roulette_wheel_selection( const std::vector< double >& probabilities, double rnd );
+
+  void
+  set_structural_plasticity_gaussian_kernel_sigma( double sigma )
+  {
+    structural_plasticity_gaussian_kernel_sigma_ = sigma;
+  }
+  /**
+   * Global list of neuron IDs used for structural plasticity computations.
+   */
+  std::vector< int > global_ids;
+
+  /**
+   * Global list of neuron positions used for spatial computations in
+   * structural plasticity.
+   */
+  std::vector< double > global_positions;
+
+  /**
+   * Standard deviation parameter for the Gaussian kernel used in
+   * spatial probability calculations.
+   */
+  double structural_plasticity_gaussian_kernel_sigma_;
+
 
 private:
   /**
@@ -195,6 +272,23 @@ class SPManager : public ManagerInterface
    */
   double structural_plasticity_update_interval_;
 
+  /**
+   * Dimentionality of the neuron positions
+   */
+  int pos_dim;
+
+  /**
+   * List of precomputed probabilities for neuron connections, indexed
+   * by neuron pair indices for efficient lookup.
+   */
+  std::vector< double > probability_list;
+
+  /**
+   * Flag indicating whether connection probabilities should be cached
+   * for performance optimization.
+   */
+  bool structural_plasticity_cache_probabilities_;
+
   /**
    * Indicates whether the Structrual Plasticity functionality is On (True) of
    * Off (False).
@@ -202,6 +296,7 @@ class SPManager : public ManagerInterface
   bool structural_plasticity_enabled_;
   std::vector< SPBuilder* > sp_conn_builders_;
 
+
   /**
    * GrowthCurve factories, indexed by growthcurvedict_ elements.
    */
@@ -228,7 +323,12 @@ SPManager::get_structural_plasticity_update_interval() const
 {
   return structural_plasticity_update_interval_;
 }
+inline double
+SPManager::get_structural_plasticity_gaussian_kernel_sigma() const
+{
+  return structural_plasticity_gaussian_kernel_sigma_;
+}
 
 } // namespace nest
 
-#endif /* #ifndef SP_MANAGER_H */
+#endif /* #ifndef SP_MANAGER_H */
\ No newline at end of file
diff --git a/pynest/nest/__init__.py b/pynest/nest/__init__.py
index 0aa91fd22b..de63103095 100644
--- a/pynest/nest/__init__.py
+++ b/pynest/nest/__init__.py
@@ -300,6 +300,20 @@ def __dir__(self):
         ),
         default=10000,
     )
+    structural_plasticity_gaussian_kernel_sigma = KernelAttribute(
+        "double",
+        (
+            "Defines the sensetivity of distance dependence in structural plasticity"
+        ),
+        default=-1,
+    )
+    structural_plasticity_cache_probabilities = KernelAttribute(
+        "bool",
+        (
+            "Controls caching of distance dependent probabiltiies in structural plasticity"
+        ),
+        default=False,
+    )
     growth_curves = KernelAttribute(
         "list[str]",
         "The list of the available structural plasticity growth curves",
diff --git a/testsuite/cpptests/run_all.cpp b/testsuite/cpptests/run_all.cpp
index fcc0868731..99b9390f2c 100644
--- a/testsuite/cpptests/run_all.cpp
+++ b/testsuite/cpptests/run_all.cpp
@@ -26,6 +26,7 @@
 
 // Includes from cpptests
 #include "test_block_vector.h"
+#include "test_distance_dependent_structural_plasticity.h"
 #include "test_enum_bitfield.h"
 #include "test_parameter.h"
 #include "test_sort.h"
diff --git a/testsuite/cpptests/test_distance_dependent_structural_plasticity.h b/testsuite/cpptests/test_distance_dependent_structural_plasticity.h
new file mode 100644
index 0000000000..d95bed9a84
--- /dev/null
+++ b/testsuite/cpptests/test_distance_dependent_structural_plasticity.h
@@ -0,0 +1,125 @@
+/*
+ *  test_distance_dependent_structural_plasticity.h
+ *
+ *  This file is part of NEST.
+ *
+ *  Copyright (C) 2004 The NEST Initiative
+ *
+ *  NEST 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 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef TEST_DISTANCE_DEPENDENT_H
+#define TEST_DISTANCE_DEPENDENT_H
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+// C++ includes:
+#include <cmath>
+#include <stdexcept>
+#include <vector>
+
+// Includes from nestkernel:
+#include "../nestkernel/random_manager.h"
+#include "../nestkernel/sp_manager.h"
+
+namespace nest
+{
+
+/**
+ * Test cases: Distance-dependent connection methods in SPManager
+ */
+BOOST_AUTO_TEST_SUITE( test_distance_dependent )
+
+BOOST_AUTO_TEST_CASE( test_gaussianKernel )
+{
+  SPManager sp_manager;
+
+  // Test for zero distance
+  std::vector< double > pos1 = { 0.0, 0.0 };
+  std::vector< double > pos2 = { 0.0, 0.0 };
+  double sigma = 1.0;
+
+  double expected = 1.0;
+  BOOST_REQUIRE_CLOSE( sp_manager.gaussian_kernel( pos1, pos2, sigma ), expected, 1e-6 );
+
+  // Test for unit distance
+  pos2 = { 1.0, 0.0 };
+  expected = std::exp( -1.0 );
+  BOOST_REQUIRE_CLOSE( sp_manager.gaussian_kernel( pos1, pos2, sigma ), expected, 1e-6 );
+
+  // Test for negative sigma (will compute as if sigma were positive)
+  sigma = -1.0;
+  double result = sp_manager.gaussian_kernel( pos1, pos2, sigma );
+  expected = std::exp( -1.0 ); // Same as sigma=1 since squared value is used
+  BOOST_REQUIRE_CLOSE( result, expected, 1e-6 );
+}
+
+BOOST_AUTO_TEST_CASE( test_get_neuron_pair_index )
+{
+  SPManager sp_manager;
+
+  // Test with valid IDs
+  BOOST_REQUIRE_EQUAL( sp_manager.get_neuron_pair_index( 1, 3 ), 3 );
+  BOOST_REQUIRE_EQUAL( sp_manager.get_neuron_pair_index( 3, 1 ), 3 );
+
+  // Test with same IDs
+  BOOST_REQUIRE_EQUAL( sp_manager.get_neuron_pair_index( 5, 5 ), 14 );
+}
+
+
+BOOST_AUTO_TEST_CASE( test_global_shuffle_spatial )
+{
+  SPManager sp_manager;
+
+  // Initialize test data
+  std::vector< size_t > pre_ids = { 1, 2 };
+  std::vector< size_t > post_ids = { 3, 4 };
+
+  sp_manager.global_ids = { 1, 2, 3, 4 };
+  sp_manager.global_positions = {
+    0.0,
+    0.0, // Neuron 1
+    1.0,
+    0.0, // Neuron 2
+    0.0,
+    1.0, // Neuron 3
+    1.0,
+    1.0 // Neuron 4
+  };
+
+  sp_manager.structural_plasticity_gaussian_kernel_sigma_ = 1.0;
+
+  std::vector< size_t > pre_ids_results;
+  std::vector< size_t > post_ids_results;
+
+  sp_manager.global_shuffle_spatial( pre_ids, post_ids, pre_ids_results, post_ids_results );
+
+  // Check result sizes
+  BOOST_REQUIRE_EQUAL( pre_ids_results.size(), 2 );
+  BOOST_REQUIRE_EQUAL( post_ids_results.size(), 2 );
+
+  // Ensure no self-connections
+  for ( size_t i = 0; i < pre_ids_results.size(); ++i )
+  {
+    BOOST_REQUIRE_NE( pre_ids_results[ i ], post_ids_results[ i ] );
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+} // namespace nest
+
+#endif /* TEST_DISTANCE_DEPENDENT_H */
diff --git a/testsuite/pytests/test_sp/test_sp_manager.py b/testsuite/pytests/test_sp/test_sp_manager.py
index 4631b613ab..799bbd410a 100644
--- a/testsuite/pytests/test_sp/test_sp_manager.py
+++ b/testsuite/pytests/test_sp/test_sp_manager.py
@@ -22,6 +22,7 @@
 import unittest
 
 import nest
+import numpy as np
 
 __author__ = "naveau"
 
@@ -48,6 +49,8 @@ def setUp(self):
             "rate_connection_delayed",
             "rate_connection_delayed_lbl",
         ]
+    HAVE_GSL = nest.ll_api.sli_func("statusdict/have_gsl ::")
+
 
     def test_register_synapses(self):
         for syn_model in nest.synapse_models:
@@ -116,7 +119,7 @@ def test_getting_kernel_status(self):
         assert sp_synapses["post_synaptic_element"] == "Den_ex"
 
         assert nest.structural_plasticity_update_interval == 10000.0
-
+    
     def test_synapse_creation(self):
         for syn_model in nest.synapse_models:
             if syn_model not in self.exclude_synapse_model:
@@ -143,6 +146,166 @@ def test_synapse_creation(self):
                 assert len(nest.GetConnections(neurons, neurons, syn_model)) == 20
                 break
 
+    
+    def test_synapse_creation_distance_dependent(self):
+        for syn_model in nest.synapse_models:
+            if syn_model not in self.exclude_synapse_model:
+                nest.ResetKernel()
+                positions = np.random.uniform([0, 0], [10, 10], size=(2, 2))
+
+                syn_dict = {"synapse_model": syn_model, "pre_synaptic_element": "SE1", "post_synaptic_element": "SE2"}
+                nest.structural_plasticity_synapses = {"syn1": syn_dict}
+                neurons = nest.Create(
+                    "iaf_psc_alpha",
+                    2,
+                    {
+                        "synaptic_elements": {
+                            "SE1": {"z": 10.0, "growth_rate": 0.0},
+                            "SE2": {"z": 10.0, "growth_rate": 0.0},
+                        }
+                    }, 
+                    positions=nest.spatial.free(pos=positions.tolist())
+                )        
+                nest.structural_plasticity_gaussian_kernel_sigma = 1.0
+                nest.structural_plasticity_cache_probabilities = True
+                nest.EnableStructuralPlasticity()
+                nest.Simulate(10.0)
+                status = nest.GetStatus(neurons, "synaptic_elements")
+                for st_neuron in status:
+                    assert st_neuron["SE1"]["z_connected"] == 10
+                    assert st_neuron["SE2"]["z_connected"] == 10
+
+                assert len(nest.GetConnections(neurons, neurons, syn_model)) == 20
+                break
+
+
+    def test_structural_plasticity_with_positions(self):
+        #Test structural plasticity when neurons have positions in multiple dimensions.
+        nest.ResetKernel()
+        nest.structural_plasticity_update_interval = 1.0
+        nest.structural_plasticity_gaussian_kernel_sigma = 1.0
+
+        syn_model = "static_synapse"
+        nest.CopyModel(syn_model, "sp_synapse")
+        nest.SetDefaults("sp_synapse", {"weight": 1.0, "delay": 1.0})
+
+        nest.structural_plasticity_synapses = {
+            "sp_syn": {
+                "synapse_model": "sp_synapse",
+                "post_synaptic_element": "Den_ex",
+                "pre_synaptic_element": "Axon_ex",
+            }
+        }
+
+        num_neurons = 30
+        positions = np.random.uniform(0, 10, (num_neurons, 3))  # 3D positions
+        neuron_params = {
+            "synaptic_elements": {
+                "Axon_ex": {"z": 1.0, "growth_rate": 1.0},
+                "Den_ex": {"z": 1.0, "growth_rate": 1.0},
+            },
+        }
+
+        neurons = nest.Create("iaf_psc_alpha", num_neurons, neuron_params, positions=nest.spatial.free(pos=positions.tolist()))
+
+        nest.EnableStructuralPlasticity()
+        nest.Simulate(10.0)
+
+        connections = nest.GetConnections(neurons, neurons, synapse_model="sp_synapse")
+        self.assertGreater(len(connections), 0, "No connections were created")
+
+        # Check that neurons closer together are more likely to be connected
+        distances = []
+        for conn in connections:
+            pre_pos = positions[conn.source-1]
+            post_pos = positions[conn.target-1]
+            distance = np.linalg.norm(np.array(pre_pos) - np.array(post_pos))
+            distances.append(distance)
+
+        avg_distance = np.mean(distances)
+        self.assertLess(avg_distance, 10.0, "Average connection distance is too large")
+    
+
+    def test_distance_dependent_without_positions(self):
+        # Test if an error is raised when distance dependency is enabled without providing positions.
+        nest.ResetKernel()
+        
+        syn_model = "static_synapse"
+        nest.CopyModel(syn_model, "sp_synapse")
+        nest.SetDefaults("sp_synapse", {"weight": 1.0, "delay": 1.0})
+
+        nest.structural_plasticity_synapses = {
+            "sp_syn": {
+                "synapse_model": "sp_synapse",
+                "post_synaptic_element": "Den_ex",
+                "pre_synaptic_element": "Axon_ex",
+            }
+        }
+
+        num_neurons = 10
+        neuron_params = {
+            "synaptic_elements": {
+                "Axon_ex": {"z": 1.0, "growth_rate": 1.0},
+                "Den_ex": {"z": 1.0, "growth_rate": 1.0},
+            },
+        }
+
+        # Create neurons without positions
+        neurons = nest.Create("iaf_psc_alpha", num_neurons, params=neuron_params)
+
+        # Enable distance-dependent structural plasticity
+        nest.structural_plasticity_gaussian_kernel_sigma = 1.0
+
+        # Expecting a NEST CppException due to missing positions
+        with self.assertRaises(nest.NESTErrors.CppException) as context:
+            nest.EnableStructuralPlasticity()
+            nest.Simulate(10.0)
+
+        # Check the error message
+        self.assertIn(
+            "Undefined positions detected for neuron ID 1. Please provide valid positions or disable distance dependency.",
+            str(context.exception)
+        )
+    def test_distance_dependent_without_positions(self):
+        nest.ResetKernel()
+        
+        syn_model = "static_synapse"
+        nest.CopyModel(syn_model, "sp_synapse")
+        nest.SetDefaults("sp_synapse", {"weight": 1.0, "delay": 1.0})
+
+        nest.structural_plasticity_synapses = {
+            "sp_syn": {
+                "synapse_model": "sp_synapse",
+                "post_synaptic_element": "Den_ex",
+                "pre_synaptic_element": "Axon_ex",
+            }
+        }
+
+        num_neurons = 10
+        neuron_params = {
+            "synaptic_elements": {
+                "Axon_ex": {"z": 1.0, "growth_rate": 1.0},
+                "Den_ex": {"z": 1.0, "growth_rate": 1.0},
+            },
+        }
+
+        # Create neurons without positions
+        neurons = nest.Create("iaf_psc_alpha", num_neurons, params=neuron_params)
+
+        # Enable distance dependency
+        nest.structural_plasticity_gaussian_kernel_sigma = 1.0
+
+        # Expecting a NEST CppException due to missing positions
+        with self.assertRaises(nest.NESTErrors.CppException) as context:
+            nest.EnableStructuralPlasticity()
+            nest.Simulate(10.0)
+
+        # Check the error message
+        self.assertIn(
+            "No neurons found. Please provide positions, or disable distance dependency",
+            str(context.exception)
+        )
+
 
 def suite():
     test_suite = unittest.makeSuite(TestStructuralPlasticityManager, "test")