diff --git a/Experiment_3/Reactive Astrocyte/Experiment_3_Reactive_astrocyte_Part2 b/Experiment_3/Reactive Astrocyte/Experiment_3_Reactive_astrocyte_Part2
new file mode 100644
index 0000000000000000000000000000000000000000..f12bc0605161bd870b615e44995ca9e67299f59c
--- /dev/null
+++ b/Experiment_3/Reactive Astrocyte/Experiment_3_Reactive_astrocyte_Part2
@@ -0,0 +1,200 @@
+"""
+
+INPUT:
+output from Part1
+mitochondria information obtained from images
+
+1- read the background mesh (bg_mesh) and level set function obtained from Part1
+2- Read the dimensionless parameters computed in Part1
+3- Read Mitochondria reaction centers point and variance obtained from microscopic image
+4- Define the reaction site coordinates
+5- Compute the Gaussian functions to locate the enzymatic activity sround the points selected in point 3 and 4
+6- Locate and define the source term of GLC close to blood vessels
+"""
+
+import numpy as np
+from dolfin import *
+from cutfem import *
+from area_ls import compute_area
+from timeit import default_timer as timer
+from Enzymes import eta_gauss, Gauss_expression_3d, uniform_coordinates_enzymes, enzymes_close2mito
+
+import time
+
+date = time.strftime("%d%m%Y")
+
+
+parameters["form_compiler"]["no-evaluate_basis_derivatives"] = False
+parameters['allow_extrapolation'] = True
+
+# Start timer
+startime_main = timer()
+
+################################
+
+# Step 1 READ FROM HDF5 File #
+
+################################
+
+f = HDF5File(mpi_comm_world(), './AD_results_'+ date + '/hdfFile', 'r')
+bg_mesh = Mesh()
+f.read(bg_mesh, "bg_mesh", False)
+W = FunctionSpace(bg_mesh, "CG", 1)
+level_set = Function(W)
+f.read(level_set, 'level_set')
+f.close()
+
+
+
+################################
+
+# Step 2 read other variables #
+
+#################################
+
+volume = np.load('./AD_results_' + date + '/Volume.npy')
+x_size, y_size, z_size = np.load('./AD_results_' + date + '/mesh_sizes.npy')
+L = np.load('./AD_results_' + date + '/L.npy')
+
+############################
+
+# The fictitious mesh
+mesh = CutFEMTools_fictitious_domain_mesh(bg_mesh,level_set,0,0);
+
+# Compute Mesh -Levelset intersection and corresponding marker
+mesh_cutter = MeshCutter(mesh,level_set)
+
+# Diminishing quadrature points
+q_degree = 3
+
+# Define new measures associated with the interior domains and facets
+dx = Measure("dx")[mesh_cutter.domain_marker()]
+dS = Measure("dS")[mesh_cutter.interior_facet_marker(0)]
+dxq = dc(0, metadata={"num_cells": 1, 'quadrature_degree': q_degree})
+dsq = dc(1, metadata={"num_cells": 1})
+
+dxc = dx(0, metadata={'quadrature_degree': q_degree}) + dxq #
+
+################################################
+
+# Step 3 Read mitochondria points and variance
+
+################################################
+
+mito = np.load('Images/AD_mito4cutfem_scaled_sigma.npy')
+mito_points = mito.T
+
+################################################
+
+# Step 4 Define the reaction site coordinates: HXK, PYRK and LDH #
+
+################################################
+
+M = 140 # number of reaction sites (chosen based on mitochondria)
+np.save('./AD_results_' + date + '/M.npy', M)
+
+# select location of HXK in a sourrounding of each mitochondria site
+coordinate_enzymes_hxk = enzymes_close2mito(level_set, mito_points, 0.03, 0.03, 0.03, x_size, y_size, z_size,M)
+
+
+# Sort the location of PYRK and LDH using uniform distribution
+coordinate_enzymes_pyrk = uniform_coordinates_enzymes(level_set, M, x_size, y_size, z_size)
+coordinate_enzymes_ldh = uniform_coordinates_enzymes(level_set, M, x_size, y_size, z_size)
+
+
+list_of_enzymes = np.ones((9,M))
+
+list_of_enzymes[0] = coordinate_enzymes_hxk[0]
+list_of_enzymes[1] = coordinate_enzymes_hxk[1]
+list_of_enzymes[2] = coordinate_enzymes_hxk[2]
+
+list_of_enzymes[3] = coordinate_enzymes_pyrk[0]
+list_of_enzymes[4] = coordinate_enzymes_pyrk[1]
+list_of_enzymes[5] = coordinate_enzymes_pyrk[2]
+
+list_of_enzymes[6] = coordinate_enzymes_ldh[0]
+list_of_enzymes[7] = coordinate_enzymes_ldh[1]
+list_of_enzymes[8] = coordinate_enzymes_ldh[2]
+
+np.save('./AD_results_' + date + '/enzymes_coordinates.npy', list_of_enzymes)
+
+################################################
+
+# Step 5 Compute the sum of the Gausssian function defined on the reaction site points #
+
+################################################
+
+# Define the dimensionless variance for the gaussians
+sigma = 1./L
+np.save('./AD_results_' + date + '/sigma.npy', sigma)
+
+# Function to define the Gaussians for mitochondria (the variance was calibrated on the image)
+def Sum_Gaussian_Mito(M, coordinate_enzymes_test, sigma):
+ Gaussian = 0
+ for i in range(M):
+ Gaussian += Gauss_expression_3d(coordinate_enzymes_test[0,i], coordinate_enzymes_test[1,i],coordinate_enzymes_test[2,i], sigma[i])
+ return(Gaussian)
+
+def Sum_Gaussian(M, coordinate_enzymes_test, sigma):
+ if M == 1:
+ Gaussian = Gauss_expression_3d(coordinate_enzymes_test[0], coordinate_enzymes_test[1],coordinate_enzymes_test[2], sigma)
+ else:
+ Gaussian = 0
+ for i in range(M):
+ Gaussian += Gauss_expression_3d(coordinate_enzymes_test[0,i], coordinate_enzymes_test[1,i],coordinate_enzymes_test[2,i], sigma)
+ return(Gaussian)
+
+Gaussian_hxk = Sum_Gaussian(M, coordinate_enzymes_hxk, sigma)
+Gaussian_pyrk = Sum_Gaussian(M, coordinate_enzymes_pyrk, sigma)
+Gaussian_ldh = Sum_Gaussian(M, coordinate_enzymes_ldh, sigma)
+Gaussian_mito = Sum_Gaussian_Mito(len(mito.T[0]), mito_points[:-1], mito_points[-1])
+
+
+# Compute the adaptive parameter for the Gaussians
+eta_hxk = eta_gauss(Gaussian_hxk, mesh, mesh_cutter, dxc)
+eta_pyrk = eta_gauss(Gaussian_pyrk, mesh, mesh_cutter, dxc)
+eta_ldh = eta_gauss(Gaussian_ldh, mesh, mesh_cutter, dxc)
+eta_mito = eta_gauss(Gaussian_mito, mesh, mesh_cutter, dxc)
+
+np.save('./AD_results_' + date + '/eta.npy', np.asarray([eta_hxk, eta_pyrk, eta_ldh, eta_mito]))
+
+##############################
+
+# Step 6 Define the source term #
+
+##############################
+
+radius_influx = 0.0307
+
+x_infl, y_infl, z_infl = 0.313278, 0.520833, 0.0302059
+x_infl2, y_infl2, z_infl2 = 0.104167, 0.0, 0.0815185
+x_infl3, y_infl3, z_infl3 = 0.364583, 0.572917, 0.0402746
+
+subdomain = Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_influx, x_0 = x_infl , y_0 = y_infl, z_0 = z_infl, degree=1)
+subdomain += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_influx, x_0 = x_infl2 , y_0 = y_infl2, z_0 = z_infl2, degree=1)
+subdomain += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_influx, x_0 = x_infl3 , y_0 = y_infl3, z_0 = z_infl3, degree=1)
+
+subdomain_volume = compute_area(bg_mesh, level_set, subdomain)
+np.save('./AD_results_' + date + '/subdomain_volume_glc.npy', subdomain_volume)
+
+##############################
+
+# Step 7 Define the LAC output #
+
+##############################
+
+radius_outflux = 0.035
+
+x_out, y_out, z_out = 0.923611, 1., 0.0203796
+x_out2, y_out2, z_out2 = 1., 0.909722, 0.0135864
+x_out3, y_out3, z_out3 = 0.506944, 0.993056, 0.0543457
+x_out4, y_out4, z_out4 =0.989583, 0.510142, 0.0402746
+
+subdomain_outflux = Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_outflux, x_0 = x_out , y_0 = y_out, z_0 = z_out, degree=1)
+subdomain_outflux += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_outflux, x_0 = x_out2 , y_0 = y_out2, z_0 = z_out2, degree=1)
+subdomain_outflux += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_outflux, x_0 = x_out3 , y_0 = y_out3, z_0 = z_out3, degree=1)
+subdomain_outflux += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? 1. : 0', r=radius_outflux, x_0 = x_out4 , y_0 = y_out4, z_0 = z_out4, degree=1)
+
+
+subdomain_outflux_volume = compute_area(bg_mesh, level_set, subdomain_outflux)
+np.save('./AD_results_' + date + '/subdomain_volume_lac.npy', subdomain_outflux_volume)
\ No newline at end of file