F[part1_, part2_, a_] := Module[{i1, j1, i2, j2, nnl},

  i1 = part1[[1]]; j1 = part1[[2]];
  i2 = part2[[1]]; j2 = part2[[2]];
  nnslist = nns[i1, j1];

  p[[i1, j1, a]]*(1 - Sum[p[[i2, j2, c]], {c, 1, N}])*Exp[Sum[-(M[[a, b]])*p[[nnslist[[neighb, 1]], nnslist[[neighb, 2]], b], {neighb, 1, 4}, {b, 1, N}]]
  ]