function [Ux, Uy, M_b, M_u, Global_Area, InsideOldMask, NatNeighborPairs] = NewNodesAdjustmentBFMyosin(new_nodes, Global_Area, Ux, Uy, M_b, M_u, PosX, PosY, InsideOldMask, BoundaryOldMask, NatNeighborPairs, NumNodesX, NumNodesY, MeshSize, Front)     

for Idx = new_nodes'
    InsideOldIdx = find(InsideOldMask); 
    [~, Id] = min((PosX(InsideOldMask)-PosX(Idx)).^2+(PosY(InsideOldMask)-PosY(Idx)).^2);
    ClosestIdx = InsideOldIdx(Id);
    NeighborOfClosestIdx = NatNeighborPairs(NatNeighborPairs(:,2)==ClosestIdx,1);
    PotentialNeighbors = [];
    for neighborIdx = NeighborOfClosestIdx'
        PotentialNeighbors = [PotentialNeighbors; NatNeighborPairs(NatNeighborPairs(:,2)==neighborIdx,1)];
    end
    PotentialNeighbors = intersect(unique([ClosestIdx; NeighborOfClosestIdx; PotentialNeighbors]), InsideOldIdx);
    
    [NatNeighbor_ij, Vertices, VoronoiCells] =...
        GetVoronoiDecomposition(Idx, PosX, PosY, InsideOldMask, PotentialNeighbors, NumNodesX, NumNodesY);
    [~, Directions, MidPoints] = GetBoundaryIndependentMetrics(NatNeighbor_ij, PosX, PosY);
    [~, Area_NewNode, ~] =...
        GetBoundaryVolumes(Idx, MeshSize, VoronoiCells, Vertices, NatNeighbor_ij, Directions, MidPoints, Front,false);
    neighbors = NatNeighbor_ij(:,1);
    InsideOldMask(Idx) = true;
    BoundaryOldMask(Idx) = true;
    NatNeighborPairs = [NatNeighborPairs; NatNeighbor_ij; NatNeighbor_ij(:, [2 1])];

    neighbors_neighbors = [];
    for neighborIdx = neighbors'
        neighbors_neighbors = [neighbors_neighbors; NatNeighborPairs(NatNeighborPairs(:,2)==neighborIdx,1)];
    end
    neighbors_neighbors = intersect(unique(neighbors_neighbors), [InsideOldIdx; Idx]);

    [NatNeighbor_ij, Vertices, VoronoiCells] =...
        GetVoronoiDecomposition(neighbors, PosX, PosY, InsideOldMask, neighbors_neighbors, NumNodesX, NumNodesY);
    [~, Directions, MidPoints] = GetBoundaryIndependentMetrics(NatNeighbor_ij, PosX, PosY);
    [~, Area_Neighbors, ~] =...
        GetBoundaryVolumes(neighbors, MeshSize, VoronoiCells, Vertices, NatNeighbor_ij, Directions, MidPoints, Front,false);
    
    Area_old_Neighbors = Global_Area(neighbors);
    Global_Area(Idx) = Area_NewNode;
    Global_Area(neighbors) = Area_Neighbors;
    
    Diff_Area = Area_old_Neighbors - Area_Neighbors;
    fi = Diff_Area/Area_NewNode;
    Ux(Idx) = sum(Ux(neighbors).*fi);
    Uy(Idx) = sum(Uy(neighbors).*fi);
    M_b(Idx) = sum(M_b(neighbors).*fi);
    M_u(Idx) = sum(M_u(neighbors).*fi);
end;