diff --git a/LoadCase/Case.py b/LoadCase/Case.py
index 52cc6222858de8852857f7ca2f1454c3029f2cd7..bc54c38ee8836e35cf566c4ad3b1249795ad3130 100644
--- a/LoadCase/Case.py
+++ b/LoadCase/Case.py
@@ -13,9 +13,9 @@ kwargs = {'x': 0, 'y': 0, 'index': 0, 'modelname': 'test', 'assembly': 0,
               'connectorE': 1500000000., 'connectorMu': 0.3, 'aMTnumber': 50,
               'aMTlength': 2, 'springStiffness': 10, 'StepName': 'Standard_Buckling',
               'CompressiveLoad': 1}
-kwargs = Interaction.create_interactions(**kwargs)
+kwargs, data = Interaction.create_interactions(**kwargs)
 Step.Step(**kwargs)
-LoadsAndBCs.create_load(**kwargs)
+LoadsAndBCs.create_load(data, **kwargs)
 LoadsAndBCs.create_bc(**kwargs)
 generate_mesh.standard_mesh(**kwargs)
 
diff --git a/LoadCase/Interaction.py b/LoadCase/Interaction.py
index a7a4a6c0acb963993030c76e610d4d842cfdbbc3..2176ec5a3487f467b2ca4765a299f5955ed59ba9 100644
--- a/LoadCase/Interaction.py
+++ b/LoadCase/Interaction.py
@@ -14,16 +14,22 @@ def create_interactions(**kwargs):
 
     :param kwargs: object
 
-    :return: (object) kwargs
+    :return: kwargs
+
+    :rtype: object
+
+    :return: data
+
+    :rtype: list
     """
-    MTdata, aMTnames, kwargs = generate_assembly(**kwargs)
+    MTdata, aMTnames, kwargs, data = generate_assembly(**kwargs)
     kwargs.update(dict(ipMTnames=MTdata['MTnames'], aMTnames=aMTnames))
 
     CoupleIpMTsToSring(**kwargs)
     CoupleIpMTsToCentrosomes(**kwargs)
     CoupleAMTs(**kwargs)
     CoupleAMTsToCentrosomes(**kwargs)
-    return kwargs
+    return kwargs, data
 
 
 def CoupleIpMTsToSring(**kwargs):
diff --git a/LoadCase/LoadsAndBCs.py b/LoadCase/LoadsAndBCs.py
index 10efd8a142dc4f7bc877387e1943f019d93e9b1b..5e72d4a7f4b5f61ee96f29d1b0f51231512e7782 100644
--- a/LoadCase/LoadsAndBCs.py
+++ b/LoadCase/LoadsAndBCs.py
@@ -2,9 +2,13 @@ from abaqus import *
 from abaqusConstants import *
 
 
-def create_load(**kwargs):
+def create_load(data, **kwargs):
     """
-    Create and apply compressing load at each growing end of ipMT
+    Create and apply compressing load at each connection between connectors and ipMTs
+
+    :param data: list containing all ipMTs and connectors of the model
+
+    :type data: list
 
     :param kwargs: model parameters
 
@@ -17,35 +21,66 @@ def create_load(**kwargs):
     # Call the current assembly
     modelname = kwargs['modelname']
     a = mdb.models[modelname].rootAssembly
-    # Select the ipMTs by names and by pole
-    ipMTnames = kwargs['ipMTnames']
     StepName = kwargs['StepName']
     LoadVal = kwargs['CompressiveLoad']
-    for i in range(len(ipMTnames)):
-        if i % 2 == 0:  # Right pole
-            v = a.instances[ipMTnames[i]].vertices
-            loc = v[-1].pointOn
-            verts = v.findAt(loc, )
-            region = a.Set(vertices=verts, name='ipMTtip-' + str(i))
-            mdb.models[modelname].ConcentratedForce(name='Load-' + str(i),
-                                                    createStepName=StepName,
-                                                    region=region,
-                                                    cf3=-LoadVal,
-                                                    distributionType=UNIFORM,
-                                                    field='',
-                                                    localCsys=None)
-        else:  # Left pole
-            v = a.instances[ipMTnames[i]].vertices
-            loc = v[-1].pointOn
-            verts = v.findAt(loc, )
-            region = a.Set(vertices=verts, name='ipMTtip-' + str(i))
-            mdb.models[modelname].ConcentratedForce(name='Load-' + str(i),
-                                                    createStepName=StepName,
-                                                    region=region,
-                                                    cf3=LoadVal,
-                                                    distributionType=UNIFORM,
-                                                    field='',
-                                                    localCsys=None)
+    # Select vertices to apply load
+    # Iterate through unique ipMT names
+    mt_names = [sublist[0] for sublist in data]
+    mt_names = set(mt_names)
+    i = 0
+    for mtname in mt_names:
+        for sublist in data:
+            if mtname == sublist[0]:
+                connector_end = sublist[3]
+                connector_name = sublist[1]
+
+                # Select points that connect ipMTs with connectors
+                ConnectorEnd = a.instances[connector_name].vertices[connector_end].pointOn
+                loadPoint = a.instances[connector_name].vertices.findAt(ConnectorEnd, )
+                loadName = 'Load-' + connector_name + '-' + str(connector_end)
+                loadRegion = a.Set(vertices=loadPoint,
+                                     name=loadName)
+                if i % 2 == 0:
+                    load = -LoadVal
+                else:
+                    load = LoadVal
+                mdb.models[modelname].ConcentratedForce(name=loadName,
+                                                        createStepName=StepName,
+                                                        region=loadRegion,
+                                                        cf3=load,
+                                                        distributionType=UNIFORM,
+                                                        field='',
+                                                        localCsys=None)
+        i += 1
+
+    # ipMTnames = kwargs['ipMTnames']
+    # StepName = kwargs['StepName']
+    # LoadVal = kwargs['CompressiveLoad']
+    # for i in range(len(ipMTnames)):
+    #     if i % 2 == 0:  # Right pole
+    #         v = a.instances[ipMTnames[i]].vertices
+    #         loc = v[-1].pointOn
+    #         verts = v.findAt(loc, )
+    #         region = a.Set(vertices=verts, name='ipMTtip-' + str(i))
+    #         mdb.models[modelname].ConcentratedForce(name='Load-' + str(i),
+    #                                                 createStepName=StepName,
+    #                                                 region=region,
+    #                                                 cf3=-LoadVal,
+    #                                                 distributionType=UNIFORM,
+    #                                                 field='',
+    #                                                 localCsys=None)
+    #     else:  # Left pole
+    #         v = a.instances[ipMTnames[i]].vertices
+    #         loc = v[-1].pointOn
+    #         verts = v.findAt(loc, )
+    #         region = a.Set(vertices=verts, name='ipMTtip-' + str(i))
+    #         mdb.models[modelname].ConcentratedForce(name='Load-' + str(i),
+    #                                                 createStepName=StepName,
+    #                                                 region=region,
+    #                                                 cf3=LoadVal,
+    #                                                 distributionType=UNIFORM,
+    #                                                 field='',
+    #                                                 localCsys=None)
 
 
 def create_bc(**kwargs):
@@ -72,3 +107,5 @@ def create_bc(**kwargs):
                                          createStepName='Initial',
                                          region=region,
                                          localCsys=None)
+
+# TODO: Apply load at the points linking connectors and ipMTs
diff --git a/SpindleAssembly/assembly_random.py b/SpindleAssembly/assembly_random.py
index b83c21bc0b87686a013b58581fad0d9b292af354..031f4de0a9ed4afe6c07fca6a5ae6c3f5d8b7b94 100644
--- a/SpindleAssembly/assembly_random.py
+++ b/SpindleAssembly/assembly_random.py
@@ -24,8 +24,8 @@ def generate_assembly(**kwargs):
     create_model_assembly(**kwargs)
     add_and_position_centrosomes(**kwargs)
     aMTnames, kwargs = add_and_position_astral_mts(**kwargs)
-    MTdata, kwargs = add_and_position_interpolar_mts_and_connectors(**kwargs)
-    return MTdata, aMTnames, kwargs
+    MTdata, kwargs, data = add_and_position_interpolar_mts_and_connectors(**kwargs)
+    return MTdata, aMTnames, kwargs, data
 
 
 def create_model_assembly(**kwargs):
@@ -67,7 +67,7 @@ def add_and_position_interpolar_mts_and_connectors(**kwargs):
     MTdata, ConnectorData, data = assign_ipMTs(**kwargs)
     attach_connectors(data, **kwargs)
     kwargs.update(dict(MTdata=MTdata, ConnectorData=ConnectorData))
-    return MTdata, kwargs
+    return MTdata, kwargs, data
 
 
 def add_and_position_astral_mts(**kwargs):
diff --git a/job.py b/job.py
index 08df4b0a6dc6ab382ca501589634649c15cb23be..ddbdc82b96eaa77716833a23b7416ae1fb4404b0 100644
--- a/job.py
+++ b/job.py
@@ -36,9 +36,9 @@ kwargs = {'x'               : 0,
           'JobName'         : 'Job-1'}
 
 ''' Call model functions '''
-kwargs = Interaction.create_interactions(**kwargs)
+kwargs, data = Interaction.create_interactions(**kwargs)
 Step.Step(**kwargs)
-LoadsAndBCs.create_load(**kwargs)
+LoadsAndBCs.create_load(data, **kwargs)
 LoadsAndBCs.create_bc(**kwargs)
 generate_mesh.standard_mesh(**kwargs)
 
@@ -52,4 +52,6 @@ mdb.Job(name=name, model=modelname, description='', type=ANALYSIS, atTime=None,
         modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='',
         scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1,
         numGPUs=0)
-mdb.jobs[name].submit(consistencyChecking=OFF)
+# mdb.jobs[name].submit(consistencyChecking=OFF)
+
+# TODO: Rerun documentation generation.